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ABSTRACT 


The  availability  of  powerful  digital  computers  has  stimulated 
new  interest  in  solving  partial  differential  equations  numerically. 
Closed  difference  methods  are  discussed  in  this  thesis;  by  closed  we 
mean  that  the  solution  is  completed  in  a  predetermined  number  of  steps. 


For  completeness ,  the  difference  operators  needed  are  defined 
and  difference  series  for  derivatives  are  derived  in  this  thesis. 

Before  employing  them  in  the  numerical  solution  of  partial  differential 
equations,  we  illustrate  their  use  in  the  simpler  context  of  ordinary 
differential  equations. 


The  analytic  properties  of  partial  differential  equations  are 
discussed  briefly  since  these  properties  help  to  determine  what  numer» 
ical  technique  should  be  used  in  a  given  problem.  The  difference  anal- 
ogue  of  a  partial  differential  equation  has,  however,  its  own  propertie 
and  difficulties  -  properties  and  difficulties  not  encountered  in  the 
analytic  solution.  An  eminent  example  of  this  type  of  difficulty  is 
the  phenomenon  of  instability:  it  is  possible  to  translate  a  well- 
posed  partial  differential  equation  problem  into  a  badly-posed  differ¬ 
ence  problem.  This  phenomenon  is  discussed  briefly  in  general  terms 
and  more  fully  in  the  case  of  the  heat  equation  and  the  equation 

— — —  =  ^  .  The  latter  is  investigated  in  detail;  the.  analytic 

5x 

solution  is  obtained  and  the  numerical  solution  undertaken. 


Matrix  methods  for  solving  linear  elliptic  systems  are 
discussed  in  detail.  These  methods  are  only  valid  if  the  boundary  is 
a  rectangle  with  sides  parallel  to  the  axes,  but  a  procedure  for 
handling  non-rectanguXar  boundaries  which  retains  the  best  qualities 
of  closed  methods  is  suggested.  Means  by  which  singularities  in 
boundary  data  are  handled  are  also  considered. 
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CHAPTER  I 

INTRODUCTION 

In  the  past  twenty  years  the  fields  of  physics,  chemistry  and 
engineering  have  called  more  and  more  on  mathematics  for  solutions  of 
differential  equations  and,  particularly  in  recent  years,  for  the 
solution  of  partial  differential  equations.  The  call  has,  to  some  extent, 
been  answered  by  the  development  of  numerical  techniques.  Since  about 
1950,  the  advent  of  electronic  digital  computers  has  vastly  increased  our 
ability  to  handle  many  such  problems  in  a  reasonable  period  of  time  [1]. 

Accompanying  the  technological  advances  in  machine  design  there 
has  been  a  flood  of  new  methods  and  techniques.  Of  these  latter,  the 
closed,  as  opposed  to  the  multiple  approximation  or  relaxation,  methods 
seem  particularly  interesting  since  they  require  little  input/output  and 
simpler  programs  and  can  be  stated  in  a  more  concise  form.  The  last 
point  is  demonstrated  in  Chapter  V  where  the  boundary  value  problem  is 
discussed  in  terms  of  the  solution  of  algebraic  equations. 

In  essence,  the  numerical  solution  of  differential  equations 
is  simple  to  describe;  the  differential  equation  and  its  accompanying 
differential  conditions  are  converted  to  a  difference  equation  and  diff¬ 
erence  conditions.  The  difference  problem  is  then  susceptible  of  a 
numerical  solution  at  discrete  points  in  the  region  of  interest. 

This  rather  general  technique  will  be  fully  illustrated  in  later 
chapters  in  reference  to  a  particular  class  of  equations,  namely  partial 
differential  equations  of  second  order  in  one  dependent  and  two  independent 
variables  having  real  coefficients.  Such  equations  can  be  classified  in 


.. 


-  2  - 

one  of  three  types  (Chapter  III)  and  each  has  a  form  of  numerical  solution 
suited  to  it.  The  details  of  the  algorithm  for  finding  the  solution,  of 
course,  depends  on  the  particular  equation  and  its  conditions.  For  example 
the  algorithm  for  solving  a  boundary  value  problem  is  quite  different  if 
the  derivative  rather  than  the  function  is  specified  on  the  boundary. 

The  analysis  of  the  errors  involved  in  the  approximation  by  a 
difference  system*  of  the  differential  system  and  of  the  errors  accrued 
while  solving  the  difference  system  is  extremely  difficult  and  far  from 
complete.  This  problem  has  been  complicated  by  the  appearance  of  distur¬ 
bances  of  a  nature  difficult  to  detect  and  connect,  but  very  destructive 
in  effect  (Chapters  II  and  IV ). 

The  validity  of  the  numerical  procedures  discussed  in  this 
thesis  and  in  the  literature  rests  on  the  existence  and  uniqueness  theorems 
of  pure  mathematics.  The  existence  problem  has  received  much  study,  with 
great  success  in  the  case  of  ordinary  differential  equations  and  limited 
success  in  partials.  In  many  cases  we  must  assume  that  the  problem  is 
well  posed  (Chapter  III)  and  in  fact  has  a  solution.  In  differential 
systems  describing  the  physical  world  this  is  generally  true.  That  the 
numerical  solution  obtained  at  discrete  points  in  space  can  be  made  to 
conform  more  and  more  closely  to  the  actual  solution,  provided  the  points 
are  taken  closer  together,  is  assumed.  This  matter  is  discussed  further 
in  Chapter  II.  Where  comparison  of  the  analytic  solution  and  the  numer¬ 
ical  solution  has  been  possible  the  results  have  been  encouraging,  but  as 
yet  little  study  of  other  than  an  empirical  nature  has  been  carried  out. 
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By  system  is  meant  equation  with  conditions  (boundary  or  initial) 
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In  what  follows,  some  of  the  closed  difference  methods  of 
solving  partial  difference  systems  in  one  dependent  and  two  independent 
variables  will  be  considered.  In  Chapter  III  the  common  difference 
operators  and  their  relationships  to  derivatives  will  be  defined.  Since, 
in  many  respects,  the  difference  methods  for  solving  ordinary  and  partial 
differential  systems  are  quite  similar,  some  consideration  will  be  given 
to  the  former  in  Chapter  II.  Chapter  III  is  devoted  to  a  brief  study  of 
the  analytic  nature  of  partial  differential  equations  in  order  to  establish 
terminology.  In  Chapter  II  there  is  a  very  brief  examination  of  parabolic 
systems.  Elliptic  equations  are  treated  in  some  detail  in  Chapter  V; 
closed,  matrix  methods  of  solution  of  elliptic  systems  are  described  and 
a  method  for  handling  irregular  boundaries  is  discussed.  The  analytic 
solution  to  the  non-linear  parabolic  equation  of  Chapter  IV  is  given  in 
Appendix  I.  It  was  hoped  that  a  comparison  of  the  analytic  and  numerical 
solutions  of  this  equation  could  have  been  made,  but  due  to  lack  of  time 
this  was  not  possible. 


Most  chapters  are  preceded  by  a  short  introduction. 
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CHAPTER  II 

DIFFERENCE  OPERATORS  AND  ORDINARY  DIFFERENTIAL  EQUATIONS 
$2.1.  Difference  Operators: 

The  definition  of  the  derivative  of  a  function,  f(x),  namely 


df  .  . 

—  =  ,  1 im_ 

dx  h  — ►  0 


f(x+h)  -  f( 

h 


*1), 


suggests  the  use  of  expansions  in  differences  of  a  function  to  represent 
derivatives . 

Suppose  we  are  given  the  sequence  of  n+1  ordinates, 


f(x0)>  f(x1)»  •••>  f(xn)  ; 


corresponding  to  the  sequence  of  n+1  abscissae. 


x05  xls  •"*  xn  * 
where,  for  h  a  fixed  interval, 

x.  =  x  +  ih,  i  =  0,  1,  2,  . . . ,  n 
i  o 

For  convenience,  write 

f (x. )  =  f . 
l  x 

The  operators,  E,  A ,\J  ,  6,  p,  D;  following  the  notation  of  Sheppard,  are 
defined  as  follows: 


(2.1) 

forward  shift, 

Ef.  = 

l 

fi+l 

i 

(2.2) 

forward  difference, 

A£i  = 

£i+l 

-  £i  - 

(2-3) 

backward  difference, 

V£i 

£i  - 

fi-l  * 

(2.1*) 

central  difference, 

5fi  = 

fi4 

-  i 

i"f 

(2.5)  averaging  operation, 


(2.6)  differentiation  operation, 


Employing  these  symbols  as  operators,  most  of  the  common  formulae 
for  interpolation,  extrapolation,  numerical  integration  and  differentiation 


can  be  derived  [2].  In  this  thesis  we  are  to  be  concerned  with  differential 
systems  and  the  derivation  of  appropriate  differential  difference  formulae 
will  be  treated  in  the  next  section.  First,  some  preliminary  results  will 
be  obtained.  In  what  follows  and  in  the  next  section  the  quantities 
E,  A,  y,  6,  |i,  D,  will  be  treated  formally.  Powers  indicate  repeated 
application  of  the  operator,  as  for  example 


From  Taylor's  theorem 


Thus  we  can  write 


(2.7) 


E  =  e 


hD 


From  (2.1),  (2.2)  and  (2.7) 


(2.8)  A  =  E  -  1 


hD 


1 


e 


Using  (2.1),  (2.3)  and  (2.7)  we  obtain 


=  6  - 


(2.9)  V=  1  -  E 


Combining  (2.l)  and  (2.4)  we  have 


and  with  (2.7)  is  derived 


(2.10)  8  =  e2hD  -  e“2hD 


=  2  sinh  JhD 


From  (2.1),  (2.5)  and  (2.7) 


(2.11) 


4  = 


J(E2  +  E“2) 


1 1  -ihD 


=  £(*2 


+  e 


=  cosh  ihD  , 


) 


Since 

cosh2  -|-hD  -  sinh2  -|hD  =  1  , 

we  have,  from  (2.10)  and  (2.1l), 

(2.12)  M=(l+iS2F  . 

Many  other  relations  of  this  kind  can  be  derived  [2],  but  these  will  be 
sufficient  for  our  purposes. 

§2.2.  Derivative-difference  Relations; 

§2.2.1.  Forward  difference  formula; 


If  we  write  (2.8) 


as 
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hD  . 
e  =  1  +  A 


and  take  the  natural  logarithm  of  both  sides  we  have 


h  D  =  In  (1  +  A)  , 


or,  for  integer  power,  m, 


(2.13)  hm  D  m  =  [in  (1  +  A)] 


m 


■  (A  -  |  A2  +  i  A3  -  . . . ) 


m 


§2.2.2.  Backward  difference  formula; 

(2.9)  can  be  written  as 


ehD  =  (1  -  V)'1 


Taking  the  log  of  both  sides  and  expanding  we  obtain 


(2. lit)  hm  D1"  =  (V+  \  V  2  +  \  V3  +  •••) 


m 


§2.2.3.  Central  difference  formulae: 

Rewrite  (2.10)  in  the  form 


hD  =2  sinh  ^  £>/2 


Divide  by  5  and  raise  to  power  m  to  obtain 


(“rX-  ( 


- 1  /  \m 


slnh~^6/2  Y 

6/2  )  ■ 


Now 


sinh 


.1  rx  _i 

x  =  /  (1  +  t2)  2  dt 

o 

rx ,,  i„2  2  k  s 

-  J  (1  -  2  +  8  ”•••) 


dt 
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after  formal  application  of  the  binomial  theorem.  Hence, 
sinh  ^x  ,12  5  4 

-  =  1  -  7  x  +  7“  x  , 

x  6  40 


and  so 


(2.15)  hm  D  m  =  6 


.in 


(i .  i. 

\  24 


2  ,  ,  \m 

+  6to  5  ' 


-)• 


In  particular,  for  m  =  2 
(2.16) 


)■ 


It  should  be  noted  that  for  m  odd  the  right  hand  side  of 
(2.15)  contains  only  odd  powers  of  8  and  for  m  even  only  even  powers. 
From  the  definition,  (2.4),  of  8  it  can  be  seen  that  odd  powers  of  5 
involve  mid- tabular  points,  i  ±  1/2,  i  +  3/2  >  ...»  which,  in  general, 
are  not  known  though  they  can  be  obtained,  with  difficulty,  by  interpol¬ 
ation.  A  more  useful  odd  derivative  central  difference  formula  is  as 
follows; 


Multiply  the  right  side  of  (2.15)  above  and  below  by  p.m  to  get 
(h  D)m  =  ( (i&)m  [p”1(l  -  §4  +  5^  8^  -  .  •  •  ]m  , 
but  from  (2.12) 


Hence,  after  multiplying  the  two  series  together  we  have 

(2.17)  hm  if  =  (u5)m  (1  -  i  62  +  6u  -  ...)"  . 


This  form  involves  ordinates  at  tabular  points  only. 
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For  a  wail  behaved  function  it  is  preferable  to  use  the  central 
difference  formulae  where  possible  since  tabular  values  on  both  sides  of 
the  point  of  application  of  the  operator  are  employed. 

$2.5.  Discussion  of  Difference  Formulae; 

To  provide  justification  for  the  above  formulae  has  proven 
extremely  difficult.  Most  arguments  in  support  of  difference  approxima¬ 
tions  are  of  an  empirical  nature  and  a  great  number  of  problems  of  vary¬ 
ing  difficulty  have  been  solved  analytically  and  numerically  to  allow 
comparison  of  solutions.  The  results  of  such  comparisons  indicate  that 
we  may  use  difference  methods  with  some  confidence. 

Two  schools  of  thought,  in  the  matter  of  difference  approxima¬ 
tions,  stand  out.  Southwell  [3]  advocates  a  crude  difference  approxima¬ 
tion  and  the  use  of  a  very  small  interval,  making  the  interval  smaller 
when  there  is  any  doubt  as  to  the  validity  of  the  results,  arguing  that 
the  discrepancy  between 

and  « 
h  dx 

may  be  made  as  small  as  we  please  provided  h  is  taken  small  enough. 
Though  this  argument  seems  reasonable  it  would  be  difficult  to  sustain 
and  in  addition  creates  many  practical  difficulties  as  we  will  see  in 
Chapter  III.  Later  work  by  Fox  [4],  [5]  has  shown  that  very  high  accuracy 
can  be  obtained  with  a  coarse  interval  if  a  higher-difference  correction 
to  the  initially  crude  difference  approximation  is  introduced.  An 
example  of  the  latter  will  be  given  later  in  this  chapter.  From  a  prac¬ 
tical  point  of  view  the  value  of  the  difference  correction  is  undoubted, 
but  again  no  rigorous  justification  has  been  given, 
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It  may  be  useful  at  this  point  to  mention  briefly  some  of  the 
theorems  of  the  theory  of  approximations  which  are  relevant  here  and  to 
explain  their  inadequacy  in  the  present  context. 

The  approximation  theorems  of  Weierstrass  are  of  fundamental 
importance  in  numerical  analysis  in  as  much  as  they  prove  that  we  can 
approximate  a  continuous  function  to  any  desired  degree  of  accuracy  by 
polynomials  or  finite  trigonometric  sums.  The  theorems,  however,  are 
of  little  practical  value  since  they  do  not  provide  an  algorithm  for 
actually  obtaining  the  polynomial  or  trigonometric  sum.  By  accepting  a 
loss  of  generality  results  of  a  more  useful  nature  can  be  obtained.  A 
proof  due  to  Bernstein  is  constructive,  but  the  Bernstein  polynomials 
converge  so  slowly  to  the  function  being  approximated  that  little  prac¬ 
tical  use  can  be  made  of  them. 

The  basis  of  the  formulae  of  §2.2  and  §2.3  was  (2.7),  namely 
f(x+h)  =  e*1  f(x) 

without  a  remainder  term,  exact  for  a  limited  class  of  functions.  In 
using  these  formulae  we,  in  most  cases,  retain  only  the  first  few  terms 
thereby  incurring  truncation  errors.  Most  analyses  of  truncation  error 
involve  some  reference  to  the  remainder  term  of  the  Taylor's  series 
(see  any  standard  text  where  interpolation  formulae  are  derived,  for 
example  ref.  6).  Hartree  [6]  shows  that  (2.7)  can  be  applied  exactly 
to  polynomials,  exponentials  (with  linear  exponents)  products  and  sums 
of  products  of  exponentials  and  polynomials. 


' 


§2.4.  Difference  Equations  as  Approximations  to  Ordinary  Differential 


Equations: 

§2.4.1.  Theoretical  Considerations: 

Although  this  thesis  is  to  be  mainly  concerned  with  the  solving 


of  partial  differential  systems  by  means  of  difference  approximations,  a 

brief  inspection  of  methods  applicable  to  ordinary  differential  systems 

* 

is  helpful  since  the  generalization  of  methods  for  O.D.E's  to  two  or 
more  independent  variables  is  not  difficult  -  formally  at  least. 

Suppose  we  are  given  the  ordinary  differential  equation  in 
the  dependent  variable  y  and  independent  variable  x 


n\ 

y  )  =  0 


(2.18)  f (x,y,y ' ,y" 


•  •  •  9 


with  appropriate  conditions.  Primes  denote  differentiation  with  respect 
to  x.  We  require  that  (2.18)  be  solvable  for  the  highest  derivative. 

It  is  required  to  find  y  =  y(x)  . 


The  theory  of  O.D.E's  yields  an  extremely  valuable  theorem 


establishing  the  existence  and  uniqueness  of  the  solution  to  the  set  of 
first  order  equations 


(2.19) 


>  •  •  •  9 


with  initial  conditions 


* 


O.D.E.  -  ordinary  differential  equation. 


C||  'f  £ 
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(2.20)  y1(xQ)  =  y°  »  y2(xc)  =  y2’  yn ( xQ )  =  y°  • 

That  (2.18)  is  equivalent  to  the  set  (2.1$))  is  seen  as  follows;  solve 
(2.18)  for  the  derivative  of  highest  order  as 


n  ,  ,  n-lN 

y  =  g(x,y,y' , . . . ,y  )  . 


Let 


yl  =  y,  yg  =  y' .  yn  =  y 


n-1 


Hence  we  have  the  set  of  equations 


yl  =  y2  ’ 


y;  =  y 


3  ’ 


r  1  —  y  > 
n-  i  n 


y'n  =  s^,yvy2,...,ya )  . 


which  is  a  special  case  of  (2.18) 


We  may  state  the  theorem  as  follows;  "Given  a  differential 
system  (2.19),  ( 2 . 20 ) ;  if  two  positive  numbers  a  and  b  can  be  found 
such  that  in  the  domain  D  defined  by  the  inequalities 

x  <  x  <  x  +  a  ; 
o  —  —  o 

y°i  -  b  £  yx  <  y°  +  b*  y2  "  b  -  y2  -  y2  +  b  •  •  •  * 

yn  -  b  *  yn  *  yn  +  b 

the  functions  f ^ ,  fg,  ...»  fft  are  continuous  and  satisfy  Lipschitz 
conditions  in  the  variables  y^  y2***»  yn  within  the  ranges,  then  a 
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positive  number  6(  <  a)  can  be  determined  such  that  in  the  interval 
(xq,  XQ+S)  the  system  possesses  one  and  only  one  solution". 

The  proof  of  this  theorem  will  not  be  given  here.  Similar 
theorems  for  boundary  conditions  can  also  be  found  [7]. 

Although  the  existence  and  uniqueness  of  the  solution  to  a 

1 

* 

given  O.D.E.  and  P.D.E.  system  may  have  been  established  by  the  general 
theorems,  obtaining  the  solution  in  analytic  form  is  another  matter.  In 
addition,  even  if  such  a  solution  can  be  found  it  is  often  extremely 
difficult  to  obtain  numerical  values  from  it  (see  for  example  Chapter  IV). 

§2.4.2.  Numerical  Considerations; 

In  general  terms,  the  numerical  solution  of  an  O.D.E.  system 
consists  of  replacing  derivatives,  where  they  appear  in  the  equation  and 
conditions,  by  a  few  terms  of  an  expression  such  as  are  given  in  §2.2 
and  solving,  if  possible,  the  difference  system  so  obtained  for  the 
numerical  values  of  the  function  at  tabular  intervals.  In  non-linear 
systems  approximations  to  powers  of  derivatives  must  be  used.  Later 
in  this  chapter  we  will  see  that  the  difference  system  should  not,  in 
general,  be  of  higher  order  than  the  differential  system.  Additional 
terms  of  the  difference  formula  may  be  included  as  a  correction  after 
a  first  approximation  to  the  solution  has  been  obtained. 

In  exemplification  of  this  statement  let  us  consider 

=  g(x,y)  , 
dx 

y(xQ)  =  a 

*  P.D.E.  -  partial  differential  equation. 


' 
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Writing 

y(xQ  +  ih)  =  y±  , 

x  +  ih  =  x,  , 
o  i 

replace  dy/dx  by  the  first  term  of  (2.I3).  The  difference  system  to 
be  solved  is 


+  gC^.y^ 


9 


yQ  =  a  ,  i  =  0,  1,  2,  . . . 


After  obtaining  the  first  approximation,  '  'y  ,  we  may  include  further 
terms  of  (2.I3)  and  solve 


(2) 


(2) 


i+1 


y.  +  g(x1,(2)yi)  +  \  A2  (l)yt  -  ±  A3  {l)y±  +  ... 


(2) 


y0  = a 


i  =  0,1,2, ... 


Generally  only  one  or  two  correction  terms  are  used. 

§2.^.3-  The  Problem  of  Stability: 

In  §2.3  ve  mentioned  one  of  the  sources  of  error  incurred  in 
the  numerical  solution  of  O.D.E.  and  P.D.E.  systems,  that  is,  the  trun¬ 
cation  of  the  difference  series.  Despite  the  difficulty  of  analysing 
such  error  we  can  place  confidence  in  difference  corrections  and  the  use 
of  small  intervals.  Inherent  in  all  numerical  work  there  are  errors  due 
to  rounding  of  figures  to  a  limited  number  of  places  and  due  to  loss  of 
significant  figures  resulting  from  arithmetic  operations.  Detailed 
analyses  of  the  latter  has  been  carried  out  for  some  of  the  well  known 
procedures  in  numerical  analysis  [8],  [9]*  Except  where  such  procedures 
are  used,  the  error  analysis  of  each  problem  depends  on  the  character¬ 
istics  of  the  problem  and  more  basic  techniques  must  be  used  [6]. 
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In  addition  to  the  above  mentioned,  a  deviation  between  the 
actual  solution  to  an  O.D.E.  or  P.D.E.  system  and  the  associated  numerical 
solution  often  occurs  which  defies  correction  and  which  exceeds  the  most 
generous  round-off “error  bounds.  This  is  the  problem  of  stability.  No 
detailed  account  of  instability  in  the  solution  of  0.0. E.  systems  will 
be  given  here.  Stability  in  P.D.E.  methods  will  be  mentioned  later 
(Chapter  III).  A  simple  example  will  serve  to  indicate  the  nature  of 
instability  in  the  solution  of  O.D.E.  systems. 

Consider  the  second  order  system 

y"  -  y  =  0  , 

(2.21) 

y(o)  =  1,  y(h)  *  e°h;  h  >  0  . 

The  analytic  solution  is 

y(x)  ~  e  X  . 

Using  the  first  term  of  equation  (2.16),  we  obtain  the  difference  system 

Vi  -  <a+h2>  -  Vi  > 

(2.22) 

~h 

y q  =  i»  y  ^  =  ®  l ,  2,  ...  * 

where  x  =  nh.  y  =  y(x  ) 
n  n  n 

Table  I  shows  the  solution  of  (2.22)  for  h  =  0.5>  0.25>  0.1. 

•nh 

The  figures  in  brackets  are  the  differences  between  e  and  the 
numerical  solution  in  units  of  the  last  figure  retained. 

Note  that  the  error  increases  alarmingly  and  rapidly  and  that 
the  solution  eventually  becomes  negative,  even  when  we  decrease  h  to 
0.1.  To  discover  the  cause  of  this  behavior  consider  the  difference 


pmiati  on  . 
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TABLE  I 

SOLUTION  OF  y(x+h)  =  (2+h2)  y(x)  -  y(x-h) 


h  = 

X 

0.5 

y(x) 

h  = 

X 

0.25 

y(x) 

h  = 

X 

0.1 

y(x) 

0 

1.00000 

0 

1.00000 

0 

1.00000 

0.5 

. 60653 

0.25 

.77880 

0.1 

.90484 

1.0 

. 36U69  (319) 

0.50 

.60627  (26) 

0.2 

.81873 

(0) 

1.5 

.21402  (915) 

0.75 

.47163  (74) 

0.3 

.74081 

(1) 

2.0 

.11685  (1848) 

1.00 

.36646  (142) 

0.4 

.67030 

(2) 

2.5 

.04889  (3320) 

1.25 

.28419  (231) 

0.5 

.60649 

w 

5-0 

1 

-.00685 

1.50 

.21968  (345) 

•  • 

•  •  •  • 

•  •  • 

•  •  •  • 

1.0 

. 36762 

(26) 

3.25 

.01293  (2585) 

•  • 

•  •  •  • 

3.50 

-.00314 

3.0 

.04638 

(361) 

•  • 

4.4 

•  •  •  • 

-.00175 

Try  a  solution  of  the  form 

y  =  c11  ,  c  a  constant 
n 

If  we  substitute  this  in  the  difference  equation  of  (2.22)  and  divide 
n— 1 

through  by  c  we  obtain 

c2  -  (2+h2)c  +1=0 


which  has  the  solutions 
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c  =  1  +  ^h2  +  h  «/l+^-h2 
=  1  +  -^h2  -  h  \A-Hj7h2 


Since 


then 


Therefore  the  general  solution  is 


.  n  _  «n 
y  =  A  c  +  Be  , 
n  ’ 

where 

c  =  >  1 

Applying  the  data  of  (2.22)  we  have 


A  +  B  =  1 

and  Ac  +  Be  ^  =  e  ^ 

Hence  B  =  1  -  A 


c  -  c 

u  The  source  of  our  trouble  is  now  apparent.  The  coefficient, 

A,  of  the  increasing  solution  to  the  difference  equation  is  not  zero 
and  despite  the  smallness  of  h  this  solution  will  disturb  and  finally 
dominate  the  solution  we  are  seeking. 

In  essence,  therefore,  instability  consists  of  the  introduction 
of  large  unwanted  solutions  to  the  difference  system.  An  instable  system 
may  be  created  in  many  ways;  by  converting  a  stable  differential  system 
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to  the  difference  system,  by  using  a  difference  system  of  higher  order 
than  the  differential  system,  and  by  trying  to  solve  an  inherently 
unstable  differential  equation. 

It  should  be  noted  that  "marching"  techniques,  which  must  be 
used  for  initial  value  problems,  (a  crude  version  of  marching  was  used 
in  the  above  example)  can  be  expected  to  exhibit  instability  since  they 
"precede"  the  data  and  "predict"  the  value  of  the  solution  at  the  next 
tabular  point.  Generally  the  numerical  methods  applicable  to  boundary 
value  problems  are  stable  provided  the  O.D.E.  system  is  stable  and  the 
order  of  the  difference  equation  does  not  exceed  that  of  the  differential 
equation. 
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CHAPTER  III 

SOME  GENERAL  ASPECTS  OF  PARTIAL  DIFFERENTIAL  EQUATIONS 
§^.1.  Introduction: 

Before  discussing  In  detail  some  numerical  procedure®  for  solving 

P.B„E.  system®,  let  us  consider  these  equations  briefly  in  a  more  general 
way  in  order  to  establish  terminology  and  concepts  which  will  b®  of  use 
in  later  chapters.  The  following  is  not  meant  to  be  complete  or  rigorous. 
Some  material  which  is  found  in  the  literature  [10,11]  is  included  here 
for  completeness  of  exposition, 

§^„2,  Definition  of  a  Partial  Differential  Equation,  Existence  and 

Uniqueness: 

A  P.D.E,  of  order  k  in  the  a  independent  variables 


*!'  X2*  °*’8  Xn 

and  one  dependent  variable 

W  "  u(xx»  x2*  *  * ‘s  Xn^ 


is  a  relation  of  the  form 


(3.1) 


b* 

dx  dx 

r  s 


) 


i , r , fe , f  ^  1,^,1,,, n , 

!t/  ■"  0,  i  o  o  a  jl) 


where  F  is  a  function  of  the  independent  variables,  the  dependent 
variable  and  the  partial  derivatives  (of  all  orders  up  to  and  including 

k)  of  the  dependent  variable. 


The  order  of  the  P.D.E.  is  the  order  of  the  highest  derivatives 


The  equation  is  linear  if  F  is  linear  in  U  and  the  partial 


derivatives  of  U  with  coefficients  depending  only  on  the  independent 
variables. 


If  F  is  linear  in  the  kth  derivatives  with  coefficients 
depending  upon  x^,  x^,  ...»  x^,  U  and  perhaps  the  partial  derivatives 
of  U  up  to  order  k  -  1,  then  the  equation  is  quasi- linear,  otherwise 
it  is  non-linear. 

The  theory  of  P.D.E  '  s  is  not  as  complete  as  that  of  O.D.E  '  s. 
There  is  no  general  theorem  of  uniqueness  and  existence  comparable  to 
that  of  §2.4.1  although  some  success  has  been  achieved  in  special  classes 
of  equations  [10],  The  theorem  of  Cauchy  and  Kowalewsky  is  one  of  the 
most  general  though  it  imposes  severe  restrictions  on  the  P.D.E.  Like 
the  theorem  of  §2.4.1  the  Cauchy-Kowalewsky  proof  [12]  applies  to  a  set 
of  P.D.E  '  s  in  the  form 

d^u1  /  du1  aV  \ 

a  k  ■  fi  {  V*2 . V  3^ . 9x  k  )  • 

1  n 

i  =  1 >  2,  . . . )  n  . 

This  is  the  normal  form.  The  f  must  be  analytic  functions  of  the  n 
independent  variables  and  the  partial  derivatives  of  the  m  dependent 
variables.  Initial  data  is  given  on  the  initial  plane  =  0  and  is 

the  form 

U  (°’*2 . xn^  =  ”l,o  (x2>x;>  •••>*!,)  > 

y-  =  <P1>t  . *„> . 


where 
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1  —  1 1  2,  •  •  •  j  in  j 
t  =  1,  2,  ...»  k-1 


and  the  functions,  <p,  must  be  analytic.  Under  these  conditions  a 
unique  solution. 


exists . 

For  the  purposes  of  this  discussion  we  will  consider  only 
linear  and  quasi-linear  second  order  equations  in  two  independent  and 
one  dependent  variable,  that  is,  equations  of  the  form 


(3.2)  aU  +2bU  +  c  U  +  d  =  0  , 

'  xx  xy  yy 

where  a,  b,  c,  d  are  real  functions  of  x,  y,  U,  U  , 
denote  partial  differentiation).  The  quantities  x,  y 
arily  denote  Cartesian  co-ordinates . 


U^  (subscripts 
do  not  neeess- 


$5.5.  Solution  Using  Initial  Data  on  a  Given  Curve; 

If  O.D.E.  theory,  when  initial  values  are  given,  we  are  often 
able  to  construct  a  solution  in  the  form  of  a  Taylor's  series  about  the 
initial  point.  For  example,  given  the  system 


y"  =  f(*,y,y’), 
y(°)  -  y0>  y’(°)  -  p  ; 

we  can  obtain  y"(o)  from  the  equation.  For  y"'(o)  we  have 


,11 1 


y  = 


df 
=  ^  + 


df  , 


df 

+ 


9 


and  so  on.  A  similar  technique  can  be  applied  to  the  partial  differential 
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initial  value  problem,  where  now  the  function  and  its  partial  derivatives 
of  first  order  are  given  along  some  initial  curve. 

Let 

U=p,  U=q,  U  =  r,  U  =  s,  U  =  t  . 
x  y  xx  xy  yy 

Then  (3. 2)  can  be  written  as 

ar  +  2bs  +  ct  +  d  =  0  . 

Let  i  be  some  parameter  on  the  curve  L  in  x,  y,  U  «*  space  on  which 
we  are  given  the  data 


(3.3)  x  =  x(i),  y  =  y(i),  U  =  U(i),  p  =  p(i),  q  =  q(i)  . 


Let  us  attempt  to  construct  a  solution  to  the  P.D.E.  system  (3.2),  (3.3) 
by  a  Taylor's  series  expansion  about  some  point  (XQ»  yQ»  Uq)  on  L: 

u(x,y)  =  Do  +  (x-xo)po+  ( y-y0)  %  +  (*-x0)  % 

+  (x°xo)(y“yo)so  +  (y=yo)2to  +  ... 

For  this  we  must  find  the  higher  order  derivatives  of  U  at  (x  ,y  ) 

00 

(we  assume  that  derivatives  of  all  orders  of  U  do,  in  fact,  exist). 


Now 


dp  =  dx  +  I?  dy 


=  r  dx  +  s  dy 


<LE  _ 

d£ 


Thus 
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Similarly 


dfl  = 
di 


s 


dx  ^  dy 

Ti  +  tTi 


These  equations,  together  with  (3.2),  give  us  three  equations  in  the 
three  unknowns  r,  s,  t;  namely: 

(3-M  Az  =  f  , 


where 


di 


0 


dx  dy 
di  di 


a 


2b 


c 


9 


r 

di 

s 

and  f  -• 

ia 

di 

t  _ 

“d 

We  can  determine  further  derivatives  of  U  by  differentiating  success¬ 
ively  (3.4).,  For  example 


dz  df  &A 

A  3“  =  3“  -  z 
ox  ox  ox 


.  dz  ^f  ^A 

A  57  =  5^  "  57  2 


from  which  we  can  determine 


dz 

57  = 


U 

U 

XXX 

xxy 

U 

xxy 

c^z 

U 

xyy 

»  57  = 

U 

u 

xyy_ 

yyyj 

Obviously,  we  can  continue  this  procedure  indefinitely  unless  we  are 
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stopped  at  the  beginning.  Suppose  that  the  determinant  of  A  is  zero 
at  the  point  on  L  that  we  are  considering,  that  is,  suppose  that 

N  -  *(£)2-2»3f  £  *«(£>  -  o  • 

This  condition  can  be  written  as 


(3-5)  a(y’)2  -  2by 5  +  c  =  0  ,  y 

which  yields  two  values  of  the  slope. 


_ 

"  dx  * 


y  dx  * 


namely 


Thus,  depending  on  whether 


A  >  0,  =  0,  <  O', 

equation  (3° 5)  defines  two  families  of  curves  in  the  x-y  plane,  one 
family  or  no  real  curves  respectively.  Notice  that  these  curves  depend 
only  on  the  coefficients,  a,  b,  c,  of  equation  (3.2).  These  are 
called  the  characteristic  curves  of  the  equation. 

Two  points  should  be  made  cle&r:  our  Taylor's  series  method 
fails  if  we  proceed  in  either  of  the  directions  defined  by  (3.5),  and 
we  cannot  obtain  a  solution  by  this  method  no  matter  what  direction 
we  go  if  the  slope  of  L  coincides  with  either  of  these  directions.  The 
characteristic  curves  (if  real)  seem  to  form  a  natural  boundary  across 
which  the  solution  may  be  discontinuous. 

It  must  not  be  supposed  that  data  at  one  point  serves  to 
determine  the  solution  throughout  the  x,y=-plane.  Since  the  character- 
istics  bound  the  region  of  validity  of  the  Taylor's  series  solution 
about  a  point  on  the  initial  curve,  all  points  on  this  curve  must  be 
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used  to  obtain  the  entire  solution. 


§5.^.  Classification  of  Partial  Differential  Equations: 

A  further  property  of  the  discriminant, 

A  =  b2  -  ac  , 


is  that  its  sign  determines  to  which  of  the  three  standard  forms 
(hyperbolic,  parabolic  or  elliptic)  equations  of  type  (3*2)  can  be 

reduced. 

Let  us  transform  (3.2),  by  a  change  of  independent  variables, 
into  a  form  in  which  the  second  order,  cross  derivative  does  not  appear, 
that  is, 


Let 


7  U65  +  6  +  d'  =  0  » 


where  7,&,d*  are  functions  of  ^,q,U,U^,U^ 


(3.6)  Cx  =  <Pb  »  \ 


£y  =  ,  T)  =  if  (^2-a)  , 


where  A^  and  A^  are  the  roots  of 


(3.7)  A^  -  (a+c)  A  -  A  =  0  , 


and  <p,  "&/  are  chosen  so  that 


^xy  ^yx  *  T''”t  ^ 


xy  yx 


du  „  du  du 

3J  -  k  5|  +  ix  3^  • 


Now 
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I 


Hence, 


d2U 


5P  “  (US  +  2^x  uc  +  (\)2  D 


ft 


=  (^)2  +  2#b2  +  fyb)2^ 


and 


d2!! 


5X  =  £  ^  U„,  +  (5  T)  +  £  T]  )  U„  +  TJ  T)  U 

oxo y  bx  *y  ££  v  by  x  bx  y'  £t)  'x  'y 


ryn 


=  cp^T^-a)  +  cp^b(>v1-2ai.A2)  +  f2b(A2>a)  U 


T]T) 


Similarly 


=  <p2(\~a)2  +  2cp^  (A1-a)(A2-a)  U.  +  ^2(A2>a)2  U 


TH 


Substitute  these  last  relations  into  (3.2). 
The  coefficient  of  is 


\i  m  2tyj  [ab2  +  b2(A1»2a+A2)  +  c(  A^a)  (A2-a)  ] 

=  2cpj (/  [eA^A  4-  (A1+A2)(b2=ac)  +  a(ac=>b2)  ]  „ 


but  from  (3.6) 


A,  A^  =  *»  A  j,  A,  +  A„  =:  a  +  c 


1  2 


12 


Hence, 


4  =  2cp^A  [“C+a+e-a] 


The  coefficient  of 


tt 


7  =  cp2  [ab2  4-  2b2  (A^-a)  4-  c(A^-a)2] 


The  coefficient  of  U  is 

T)T\ 


. 
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5  =  ^r2  [ab2  +  2b2(A2>a)  +  c(A2~a)2]  . 

Ignoring  for  the  moment  the  distinction  between  A^  and  »  we  can 
write  7  /  cp2  and  5  /  if/2  as 

€  =  ab2  +  2b2  (A~a)  +  c(A-a)2  . 

From  (3.7)  we  have 

b2  -  (A-a)(A-c)  . 

Replace  b  in  the  first  term  of  the  expression  for  e  by  the  above 
expression  and  we  obtain 

e  =  a(A=a)(A“c)  +  2b2(A-a)  +  c(A-a)2 

=  (A=a)  [A(a+c)  +  2(b2-ac)] 

=  (A-a)  [A(A1+A2)  -  2  A^]  - 

Thus 

7  /  <PS  =  (*j-a)  \  > 

8  /  fa  =  -(A2-a)  (Ax-A2) 

and  so 

- MM  \  \  ■ 

qp^2 

Now 

(Ara)(A2-a)  =  Wa(W  +  ^ 

p 

=  "b2  +  ac  -  a(a+c)  +  a 
=  -b2  . 
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Therefore 


t 

where  we  choose 


•v  a+c 
\ - 


rSfa  =  -b2  (\  -  \)s  A  , 

<p  H/  12 


■_f(a+c)g  +  4A^ 


a  =  a+_c_+JIa+c)2  +  UaF 


There  are  three  cases  to  consider. 
Case  (i);  A  =  0  . 


Hence, 

\  =  °»  \  ~  a+c 

and  so 


7  =  0, 

S  =  f2  c(a+c)£  . 

Case  ( ii) :  A  >  0  . 

Signs  of  7  and  6  are  opposite. 

Case  (iii):  A  <  0  . 

Signs  of  7  and  5  are  the  same. 

Therefore,  if  we  write 

a2  =  „2  | (Y«)  \  (Aj-y  I  • 

e2 .  KVa)  \  <W  I 

equation  (3.2)  can  be  reduced  to  one  of  the  following  standard  forms 


1 2  (  \2  d2''-1  , 

^  c(a+c)^  ^ps  +  d' 


2  ^2U  (j2  ^2U  ,  , 

■  _  “  P  +  d 


a 


2  ^2U  n  2 

a  5^  +  P  S^i  +  d 


o  ,  a  =  o  ; 

o  ,  a  >  o  ; 

o  ,  a  <  o  . 
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By  comparing  these  equations  with  the  algebraic  equations 

y2  =  const .  x  , 
a2*  -  p^2  =  const,, 

a2x2  +  p^2  =  const,  , 

we  name  the  standard  forms  parabolic,  hyperbolic  and  elliptic  respectively. 
Notice  that  parabolic  P.D.E  '  s  have  one  real  characteristic  curve,  hyper¬ 
bolic  equations  have  two  and  elliptic  equations  have  none. 

As  an  example  consider  the  equation 

U  +  y  U  =  0  . 

xx  yy 

Now 

A  =  b2  -  ac  =  =y 

Thus,  the  P.D.E,  is  hyperbolic  in  the  half  plane 

y  <  o  , 

elliptic  in  the  half  plane 

y  >  0  , 

and  parabolic  along  the  line 


y  =  o  . 


The  characteristic  curves  are  given  by 


( 


& 

dx 


+ 


y  =  o  , 


which  integrates  to 


$  \ 
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y  =  -  ^  (x  +  e)2  ,  e  being  a  constant. 

These  curves  are  real  only  if  y  <  0  . 

It  does  not  follow  that  transforming  a  given  P.D.E.  to  standard 
form  provides  the  best  method  for  practical  solution;  indeed,  such  a 
transformation  may  unnecessarily  complicate  the  problem  since  the  given 
data  must  also  be  transferred.  The  practical  difficulty  of  solving 
equations  (3 .6)  for  £(x,y)  and  q(x,y)  and  equation  (3.5)  for  the 
characteristics  should  not  be  ignored,  since  if  the  equation  is  not  linear 
this  involves  finding  the  solution  to  the  equation  as  well.  Several 
methods  for  solving  quasi- linear  systems  in  terms  of  characteristics  have 
been  developed  [10]  which  involve  the  simultaneous  determination  of 
t) ,  U,  and  .  Corresponding  numerical  procedures  exist  [13]. 

S3. 5.  Data  for  Partial  Differential  Equations; 

In  the  preceding  sections  of  this  chapter  we  found  that 
hyperbolic,  parabolic  and  elliptic  P.D.E  's  differed  from  each  other  in 
the  possession  of  characteristic  curves.  It  is  reasonable  to  ask  if 
there  are  other  properties  which  distinguish  these  equations  from  each 
other . 

In  a  practical  sense  it  is  quite  meaningless  to  talk  of  solving 
a  P.D.E.  without  making  reference  to  certain  auxiliary  restrictions  on 
the  dependent  variable  and,  perhaps,  its  derivatives.  We  will  define 
boundary  conditions  to  mean  data  given  on  a  closed  curve  in  the  plane 

of  the  dependent  variable,  and  initial  conditions  to  mean  data  given  on 

•  « 

an  open  curve.  We  might  ask  if  each  type  of  P.D.E.  requires  boundary  or 
initial  data  in  order  that  the  problem  be  well-posed  in  the  sense  that 
the  solution  should  exist,  be  uniquely  determined  and  depend  continuously 
on  the  data. 
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Let  us  consider  a  few  examples  from  physics,  assuming  that 
P.D.E.  systems  resulting  from  the  correct  application  of  physical  laws 
represent  well-posed  problems. 


The  wave  equation  in  one  space  dimension; 

The  equation  of  motion  of  a  string  of  finite  length,  under 
uniform  tension  T,  with  uniform  density  e,  undergoing  small 
amplitude  motion  is 


afu 

dx2 


1  dfU 
c2  at2 


=  T  /  e  >  0  , 


where  x  is  the  space  co-ordinate  and  t  is  the  time  co-ordinate. 
Clearly,  the  equation  is  hyperbolic. 


Auxiliary  conditions  might  be  obtained  by  fixing  the  ends  of 
the  string  and  displacing  it  from  the  point  of  rest  (x^o)  into  some 
configuration,  say  f(x),  then  releasing  it.  If  the  ends  of  the  string 
are  the  points 


x  =  a,  x  =  b 

the  auxiliary  conditions  can  be  stated  as  follows? 

U(a,t)  =  U(b9t)  =  0  , 

U(x,0)  =  f(x)  , 

Ut(x,0)  =  0  , 
t  >  0  ,  a<x<b. 

Conditions  such  as  these  are  initial  since  they  are  given  on  an  open 
curve  defined  by 

x  =  a,  x  =  b, 


t  =  0  . 


■ 

. 
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Hence,  we  might  suspect  that  hyperbolic  equations  require  initial  conditions. 
§3»5<>2.  Poisson's  equation; 

Suppose  we  are  given  a  long  rectangular  box  of  width  a, 
height  b,  filled  with  a  dielectric  medium.  The  opposing  sides  of  the 

box  are  at  potential  ¥. ,  V  ,  V  ,  V.  respectively  and  the  charge  density 

12  3  ^ 

in  the  medium  is  f(x,y)„  If  U  is  the  potential  in  the  box,  then 
D  =  U(x,j) 


with  no  dependence  on  z  .  The  equation  for  the  potential  is 

3%  d2U  ,  . 

“  +  “ “  =  f(x,y)  . 

dx2  dy2 


If  the  box  is  symmetric  about  the  origin,  the  auxiliary  conditions  are 


u(f,  7)  -  Vx  ,  U(-|,  y)  =  V2  , 
U(x,  |)  =  V  ,  U(x,  -|)  =  . 


The  equation  is  elliptic  and  the  data  is  of  boundary  type.  We  anticipate 
that  well-posed  elliptic  problems  require  boundary  data. 


§3° 5° 3°  The  heat  equation; 

Let  U(x,t)  be  the  temperature  at  time  t  and  at  point  x 
of  a  narrow,  uniform  rod  lying  along  the  x  °  axis.  If  there  are  no  heat 
sources  or  sinks  in  the  rod,  the  equation  governing  the  temperature  dis¬ 
tribution  is 


^x2 


t  >  0  ,  x  >  0  , 


where  a2  is  the  constant  of  heat  diffusivity  of  the  rod.  This  equation 


is  parabolic. 


' 
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Suppose  the  rod  is  infinitely  long  and  initially  has  the  temp¬ 
erature  distribution  f(x).  Then  the  appropriate  condition  is 

U(x,0)  =  f(x) 

which  is  an  initial  condition.  Hence,  parabolic  systems  appear  to  be 
initial. 


In  summary  then,  we  are  led  to  suspect,  by  considering  physical 
problems,  that  we  1 Imposed  hyperbolic,  elliptic  and  parabolic  systems  are 
initial,  boundary  and  initial  value  problems  respectively.  Examples  can 
be  given  [12]  to  show  that  applying  initial  data  to  elliptic  equations 
or  boundary  data  to  hyperbolic  or  parabolic  equations  will  result  in 
systems  which  are  not  well-posed. 

§5.6.  Remarks  on  the  Numerical  Solution  of  Partial  Differential  Systems: 
§3.6.1.  Difference  approximations  to  parjbial  differential  systems; 

Although  the  numerical  solution  of  P.D.  systems  has  received 
much  attention  since  early  in  the  1940”s  [1,3, li+]  (see  esp.  [1]  for  a 
comprehensive  bibliography)  little  of  a  general  nature  can  b®  said. 

In  this  section  a  procedure  for  reducing  equations  of  type 
(3.2)  to  difference  equations  will  be  discussed.  Essentially,  this 
reduction  is  the  same  as  that  of  an  O.D.  system  except  that  the  solution 
is  sought  at  discrete  points  in  a  two-dimensional  region  rather  than  on 
a  line. 


For  the  moment,  assume  that  the  data  curve  is  a  rectangle  com¬ 
posed  of  the  axes  and  the  lines  x  =  1,  y  s  1  if  it  is  a  boundary  prob¬ 
lem.  If  the  problem  is  of  initial  type,  assume  the  initial  curve  is  an 
open  rectangle  bounded  on  three  sides  by  the  axes  and  the  line  x  =  1. 


>  ?. 


■=  3k  " 

By  a  linear  change  of  variables  we  can  always  obtain  this  configuration 

provided  the  data  curve  is  rectangular.  Non-rectangular  curves  present 

[ 

a  difficulty  to  be  discussed  later.  Let  us  also  assume  that  we  are  to 
find  solutions  to  the  initial  value  problems  in  the  range  0  <  y  <  1. 

Let  us  examine  in  detail  the  procedure  for  setting  up  and 
solving  the  difference  system  corresponding  to  the  P.D.  system 


aU  +  bU  +  cU  +  d  =  0  , 
xx  xy  yy 

0  <  x  <  1  ,  0  <  y  <  1  . 

If  the  system  is  initial,  let  the  data  curve  be 
x  =  0  ,  x  =  1  ,  y  =  0  . 

If  the  system  is  boundary,  consider  the  data  curve 


x  =  0  ,  x  =  1  ,  y  -  0  3  y=l„ 

A  mesh  or  grid  is  formed  in  the  region  of  solution  composed  of 
vertical  and  horizontal  lines  through  the  points 


where 


h 

x 


and 


i  =  0,  1,  2,  . . . ,  m  , 


j  —  0,  1,  2, 


o  &  o  5 


in  o 


. 


and  let  corresponding  expressions  for  a,  bs  c,  d  be  denoted  a^  , 
bijs  Cij’  dij  * 


Partial  derivatives  with  respect  to  x  and  y  are  replaced 
by  difference  formulae  drawn  from  §2.2  with  differences  taken  in  the  x 
and  y  directions  respectively.  Subscripts  placed  on  difference  oper¬ 
ators  will  denote  the  direction  in  which  differencing  is  to  occur.  The 
difference  equation  can  now  be  written  down  at  each  point  of  the  region 
of  solution,  as  for  example, 


a  5  U. .  +  b 
ij  x  ij  ij 


4  S  ( p  5  U.  . )  +  e .  .  5 
x  x  y  y  ij  ij  1 


ij 


+  d. .  =  0 
ij 


If  the  equation  is  not  linear  we  can  say  little  more  in 
general  since  completion  of  the  algorithm  depends  on  the  particular 
equation  being  considered.  It  is  common  practice  to  linearize  the 
difference  equation  by  the  device  of  replacing  non-linear  terms  by  linear 
approximations.  A  crude  first  approximation  to  the  solution  can  be 
obtained  in  this  way.  After  the  first  solution  is  obtained  we  can 
linearize  the  equation  by  using  this  estimate  in  much  the  same  manner 
as  difference  corrections  are  applied.  For  example,  a  term  such  as 


.  might  be  replaced  fey  .  ,  where 


1) 


ij 


U.  .  is  the 
ij 


2N 


known  first  estimate  and  '‘"‘'U,  .  is  to  be  found.  Difference  correct- 

ij 

ions  can  be  applied  as  well.  An  example  of  the  above  procedure  is  given 
in  Chapter  IV  where  the  e« 


Sx2 


og 

St 


is  considered. 
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Even  if  the  P.D.E.  is  linear  a  complete  discussion  of  the 
algorithms  for  solving  the  P.D.  system  would  necessitate  considering 
all  possible  data  conditionss  each  requiring  the  difference  equation  in 
a  different  form,  A  particular  example  is  given  below. 

§3.6.2.  Example; 

To  illustrate  what  may  be  called  the  general  procedure  for 
setting  up  the  closed  difference  solution  to  a  P.D.  system,  let  us 
consider  the  following  example; 


(3.8) 

Given  the  linear s  hyperbolic  P.D.E. 

aU  +  bU  +  cU  4-  dU  4-  eU  +  fU  4-  g  =  0 
xx  xy  yy  x  y 

and  the  data 


(3.8a) 

U(0,y)  =  A(y)  , 

(3.8b) 

U(l,y)  =  B(y)  , 

(3.8c) 

U(xsO)  --  C(x)  , 

(3„8d) 

(  f  )y.O  -  D(x)  ’ 

0<x<ls  0  <  y  <  1  s 

where  a,  bs  cs  d9  es  f9  gs  A,  Bs  Cs  D  are  well-behaved  functions. 
Assume  that  the  data  has  no  discontinuities  at  the  corners. 

We  require  that  the  equation  be  hyperbolic  throughout  the 
region  of  solution  and  that  the  line  y  =  0  not  be  a  characteristic. 
It  is  assumed  that  the  solution  exists „  is  unique  and  that  all  deriva¬ 
tives  of  U  with  respect  to  x  and  y  exist. 


Using  (3.8)s  (3.8c)  and  (3«8d)  we  can  obtain  the  first  row 
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solution,  U  ^  ,  by  a  Taylor 9  s  series^ 

Uil  “  U(x’hy>  "  U(x-°>  +  hy  (|  )y=0+  hy  (  P  )  y=0  + 


Write  (2.16)  in  the  form 


,2n2  R2 

h  D  =8  +  c. 


and  (2.I7)  as 


hD  =  p&  +  c. 


-  i  (E-E"1)  +  c2 


From  the  latter  we  have 


h  h  U  =(46+0  _)(p  5  U+c  _U) 
x  y  xy  x  x  x2  y  y  y2 


=  (j’S  (p  8  U)  +  c  U 
x  x  y  y  5 


The  expressions  represented  by  c^s  c^s  c^  are  dif  ference<=correction 

terms „ 

2  2 

If  we  multiply  (3„8)  by  h  h  and  replace  derivatives  by 

x  y 

difference  formulae,  we  obtain 

ah2(8213+c  ,U)  +  bh  h  [p  8  (p  8  U)+c,U]  +  ch2(82B+c  ,U) 
y  x  xl  '  x  y  x  x  y  y  3  x  y  yl  ' 

+  dh  h2(p  5  U+c  U)  +  eh^  (p  8  U+c  0U)  +  fh2h2U+gh2h2  =  0 
x  y'  x  x  x2  x  y  'y  y  y2  x  y  x  y 

2 

Dividing  through  by  h^,  taking  all  correction  terms  to  the  right  and 
writing 


U(x,y)  =  U 


ij 


gives  us  the  difference  equation 


' 
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where 


a.  k282U  .  +  b  .k|i  6  (|i  5  U  )  +  c  .&2U,  .  +  d  h  k(i  5  U.  . 
ij  x  ij  ij  ^x  x'  y  y  ij'  ij  y  ij  ij  y  Kx  x  ij 

+  e  h  (i  5  U  +  f .  .h2U  .  +  g.  ,h2  =  c  U.  .  , 
ij  y  y  y  ij  ij  y  ij  ij  y  5  ij 


i  =  1,  2,  . . . ,  n  -  1  , 


J  ~  1 )  2,  •  •  •  ,  in  °  1  j 


k  = 


=  h/'x 


The  given  data  is 


U  =  A .  ,  U  .  =  B  „  j  U  =  C 
oj  J  nj  J  io  i 


and  U  ^  has  been  calculated, 


Write  the  difference  equation  explicitly  in  terms  of  the 

unknown  row,  U#  .  ,  (J=0, 1,2, . . . ,n) ,  as  follows? 

^ ,  j+i 


=  c  U 

5  ij 


-  |"a.  .k2B2' 

|_  ij  x 


V 

U.  .  +  ■ — ■* 


U 


u, 


ij  4  V  i-l»j“l  i+1, j-1 


s  \  d  „  .h  k  e  „  .h 

+  c,  / U.  .  ,  “2U  ,  )  +  XX™  |i  5  U,  .  -  U.  .  . 

ijV  i,j-l  ij J  2  x  ij  2  i,J-l 


+  h  (  f .  .U. .  +  g  . 
y  V  ij  ij  ij. 


To  simplify  notation,  write  this  as 


a. .U.  -  .  ,  +  3 .  .U.  .  .  -  a. .U.  ,  .  ,  =  V  +  c  U  . 
ij  l-l, j+1  Hij  i,j+l  ij  i+1 s J+l  i  5  ij 


Hence,  for  i  =  1,2, . . . ,n=l,  we  have  (n~l)  equations  in  the 
(n~l)  unknowns 


-  39  - 


U.  •  1 
1,  J+1 

cs 

CM 

P 

* 

j+1  ’ 

9  9  O  f 

u  .  ,  .  . 

n-1, j+1 

These 

equations , 

together  with 

the 

data,  can  be 

written  in  matrix 

form, 

(3«9) 

FU  = 

R  , 

where 

’alj 

0 

0 

0  0 

a2j 

“a2j 

0 

0  0 

F  = 

0 

• 

a3J 

9 

“CO  —  •  9  0  9 

9  9 

0  0 

9  0 

9 

0 

0 

0 

0 

a  1  .  B  ,  , 

n-1 » J  Mn-l,j_ 

. 

__ 

u  = 


u 


i,j+i 


u 


2, j+1 


U 


n-1, j+1 


R  = 


V1  +  c5  Ulj  ‘  alj  V 


v2  +  C5  U2j 


Vn-2  +  C5  Un-2, j 

Vn=l  +  C5  Un-l,j  +  an=lsj  Bj+1 


The  procedure  is  to  let 


°5  Uij  '  ° 


■ 


-  Uo  - 


and  solve  the  matrix  equation,  (3.9) >  successively  for 


J  —  1,2, .. . ,m- 1  . 


This  gives  the  first  approximate  solution, 
is  then  evaluated,  included  in  the  vector 
is  solved,  as  before,  for 


The  correction  term,  c^U^., 

5  ij 

R  and  the  corrected  equation 


j  =  1, . . . ,m-l  . 


This  procedure  is  repeated  until  successive  solutions  agree  to  the 
required  accuracy.  An  additional  condition,  usually  applied  when  P.D. 
systems  are  to  be  solved  by  difference  methods,  is  to  require  that 
differences  (of  the  unknown  function)  of  some  order  be  less  than  a 
prescribed  constant.  This  serves  to  determine  the  mesh  sizes  h  and 

X 

h  . 

y 

§3.6.3.  The  state  of  the  art: 

Though  the  discussion  of  the  preceding  two  sections  was  rather 
restricted  in  scope,  it  did  serve  to  illustrate  the  manner  in  which  P.D. 
systems  are  tackled  numerically.  We  see  that  the  problem  is  more  diff¬ 
icult  than  a  step  from  one  to  two  dimensions  would  indicate  and  it  is 
clear  that  the  state  of  the  art  is  nowhere  as  advanced  as  is  the  case  in 
O.D.E  's.  The  analysis  of  truncation  and  round-off  errors  and  stability 
is  extremely  difficult,  particularly  for  quasi  and  non-linear  systems. 

The  advent  of  digital  computers  has  made  the  solution  of 
two  dimensional  systems  common  and  three  dimensional  systems  possible, 
but  computers  have  not  cured  all  our  ills.  A  great  deal  of  time  must 
be  spent  analysing  and  preparing  the  problem  for  the  machine.  Programs 
for  solving  P.D.  systems  are  generally  very  complicated,  involving 
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either  a  great  deal  of  internal  decision-making  or  input  and  output  to 
permit  human  intervention.  The  great  strength  of  digital  computers  is 
the  ability  to  perform  arithmetic  operations  rapidly  and  it  is  only  with 
the  greatest  of  difficulty  that  a  computer  can  be  programmed  to  even 
approximate  human  decision-making  ability.  In  order  to  best  make  use 
of  human  intelligence  as  a  machine  saving  device,  considerable  communi¬ 
cation  between  the  machine  and  its  operator  is  required  and  as  yet  the 
speed  of  input/output  devices  lags  far  behind  the  operating  speed  of 
computers . 


Writing  the  difference  equation  in  matrix  form,  as  we  have 
done  in  §3*6.2,  serves  to  emphasize  the  main  drawback  to  decreasing  the 
mesh  size  to  ensure  accuracy.  The  number  of  equations  to  be  solved 
varies  as  l/h^h^  and  the  time  taken  to  solve  the  systems  of  linear 
equations  required  to  obtain  the  complete  solution  increases  as  the 
cube  of  the  number  of  equations,  that  is,  as  (l/h  h  .  Furthermore, 
the  more  equations  there  are  to  be  solved  the  greater  the  number  of 
computations  which  must  be  performed?  this  increases  round-off  error 
and  hence  more  significant  figures  must  be  carried,  which  serves  to 
slow  down  computation.  Finally,  big  matrices  of  numbers  with  many 
significant  digits  require  large  blocks  of  storage  in  the  computer. 

If  sufficient  high  speed  storage  (main  memory)  is  not  available,  slow 
access  storage  (magnetic  tape,  for  example)  and  more  complicated  methods 
of  solving  linear  equations  must  be  used.  This  again  is  costly  in  time. 


Difference-correction  techniques  permit  larger  intervals,  but 
they  also  suffer  from  certain  drawbacks.  In  order  to  calculate  differ¬ 
ences  near  the  data  curve,  values  of  the  unknown  function  outside  the 
region  of  solution  must  be  known.  These  values  are  obtained  either  by 


■ 


-  k2  - 


extrapolating  the  known  solution  or  by  solving  the  difference  equation 
with  already  computed  values  as  data.  Both  these  procedures  may  be 
risky,  especially  if  there  is  a  possibility  of  crossing  a  characteristic 
curve.  Since  we  are  to  eventually  take  differences  of  the  exterior 
values,  error  would  tend  to  be  minimized,  provided  we  know  the  solution 
to  be  continuous  across  the  data  curve.  This  is  not  always  the  case  in 
initial  value  problems  though  it  is  usually  true  in  elliptic  systems. 

The  evaluation  of  the  dif ference-correction  term  requires  a  fairly 
complicated  program  and  in  addition  memory  requirements  for  the  storage 
of  the  program  and  tables  of  differences  are  large,  hence,  auxiliary 
storage  is  usually  necessary.  All  differences  greater  than  some  pre- 
determined  value  must  be  included  in  the  correction  term.  Since  the 
magnitude  of  the  difference  depends  on  the  mesh  size,  it  may  happen 
that  the  calculation  of  the  correction  term  becomes  tedious  indeed  if 
the  mesh  size  is  too  large.  A  combination  of  decreasing- interval  and 
difference-correction  techniques  appears  to  be  the  best  approach;  for 
example,  we  could  require  that  differences  greater  than  some  given  order 
are  to  be  negligible  and  decrease  the  mesh  size  accordingly.  In  any 
case,  runs  of  two  or  three  days  of  continuous  computing  are  not  unusual. 

Non-rectangular  data  curves  make  solving  the  difference 
equation  extremely  tedious  since  the  curve  will  not  pass  through  mesh 
points.  If  the  mesh  is  adjusted  to  do  so,  then  all  grid  rectangles  are 
not  the  same  size  and  closed  difference  methods  cannot  be  applied  so 
readily.  One  method  of  dealing  with  non-rectangular  data  curves  is 
discussed  in  Chapter  V.  Sometimes  it  is  possible  to  obtain  a  rectang¬ 
ular  data  curve  by  a  change  of  dependent  variables.  In  most  cases  the 
solution  must  be  interpolated  to  the  boundary. 
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It  is  clear  that  much  basic  work  must  be  done  to  develop  more 
efficient  methods  for  solving  P.D.  systems  and  faster  computers  with 
increased  memory  capacity,  decision-making  ability  and  improved  input/ 
output  facilities  to  implement  these  methods,  before  consideration  of 
P.D.E  '  s  can  be  brought  to  any  degree  of  completion. 


I 


. 
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CHAPTER  IV 


SOLUTION  OF  THE  EQUATION 


d2^ 

dx2 


du 

ST 


$4.1.  Introduction; 

As  an,  example  of  the  application  of  difference  methods  to 
quasi- linear  equations  Richtmyer  [14]  considered  the  equation 


(t.l) 


o2]]' 

dx2 


du 

St  ’ 


t  >  0  . 


By  requiring  U  to  be  a  running-wave  solution  of  (4.1),  that  is. 


U  =  U(x  -  vt)  , 


Richtmyer  obtained  the  following  solution  in  implicit  form; 

f  (U-U  U  (U-U  )3  +  15  U2(U-U  )2  +  20  U5(U~U  ) 

4  o'  3  ov  o'  o'  o'  o'  o' 

+  5  Uq  in(U“UQ)  =  v(vt-x+x  )  , 

where  U  ,  x  and  v  are  constants.  Numerical  values  were  assigned  to 

o '  o 

these  constants,  the  equation  was  solved  for  U  by  Newton's  method  and 
compared  with  values  obtained  by  the  numerical  solution  of  (4.1).  In 
general,  the  two  solutions  agreed  closely. 

A  more  direct  solution  to  (4.1)  is  given  in  this  chapter  and 
Appendix  I.  It  was  hoped  that  a  comparison  of  analytic  and  numerical 
solutions  could  be  made  here,  but  the  task  is  difficult  and  time  consuming 
and  complete  numerical  results  are  not  yet  available.  Some  partial  results 
are  recorded  briefly  at  the  end  of  this  chapter. 


In  this  next  section  a  brief  discussion  of  the  numerical 
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properties  of  parabolic  systems  will  be  given  with  particular  reference 
to  (4.1)  and  the  heat  equation 


(4.2) 


2 

a 


dfu 

6x2 


du 

ST 


§4._2._ The  Parabolic  Difference  Equation: 

§4.2.1.  Implicit  and  explicit  forms  of  the  difference  equation; 

In  approximating  any  P.D.E.  by  a  difference  equation,  we  have 
to  decide  which  of  the  formulae  of  §2.2  to  use.  By  considering  the 
linear  heat  flow  equation,  (4.2),  we  will  see  that,  in  parabolic  systems 
at  least,  this  choice  is  crucial.  As  physical  intuition  might  tell  us, 
it  is  the  essential  difference  between  time-like  and  space-like  co-or¬ 
dinates  which  is  the  determining  factor. 


Established  practice  suggests  that  we  replace  the  second  order 


derivative  by  the  central  difference  formula,  (2.16).  The  time  deri- 
dU 


:ive,  ,  can  be  replaced  by  one  of  three  formulae:  forward  (2.13)S 

backward  (2.14),  or  central  (2.I7);  each  of  which  results  in  a  slightly 
different  difference  equation.  To  illustrate  these  three  equations 
consider  a  typical  mesh  point  in  the  region  of  solution  and  label  it 
0;  surrounding  mesh  points  are  labelled  1,  2,  3»  •>  °  °  according  to  the 
figure  below. 


(0 

1 

\S 

3 

0 

i 

( 

0— — 0 

7 

x-di 

Let 

rection  is  t< 

the  value  0 

f 

d  the  right 

f  U  at  the 

T 

and  the  positive  t 

mesh  point  i  by 

mesh  points  involved  in  the  three  difference  equations  which  follow 
appear  in  diagrams  to  the  right  of  the  equations.  Lines  joining  points 


-  46  - 


indicate  differences  involving  these  points.  Dots,  (*),  represent  knbwn 
points  and  crosses,  (x),  are  unknown  points.  Only  the  first  terms  of  the 
respective  difference  formulae  are  used. 

Using  (2.I5)  and  (2.16)  at  0  gives 


(If. 3)  ht  «2(U1  +  D?  -  2U0)  =  h2(U£  -  Uo) 


From  (2.14)  and  (2.16)  applied  at  2  we  have 
(k.h)  ht  a2(U5  +  U6  -  2Ug)  .  h2(U2  -  Bo)  . 


0 


0 


Finally,  if  (2.I7)  and  (2.16)  are  applied  at  09  we  obtain 


,5)  ht  a2(U1  +  U?  -  2Uo)  -  h|(U2  -  OjP 


3 


2 

x 


0 


3- 


Suppose  that  we  are  given  the  values  of  U  along  t  =  0  and  that  we 
are  carrying  the  solution  forward  in  the  direction  of  increasing  t, 
one  row  at  a  time.  Equation  (4.5)  is  not  suitable  in  this  case  since 
it  requires  two  rows  of  data  in  order  to  get  started.  Equation  (4. 5) 
can  be  solved  explicitly  for  U2  provided  U„,  Uq  and  ^  are  known. 
Assuming  values  at  all  the  points  . . . ,3»0,1, . . .  are  known  we  can  u3e 
(4.3)  repeatedly  to  obtain  the  next  row,  . . . ,6,2, 5» ° ° «  °  Because  the 
values  on  one  row  are  given  explicitly  point  by  point  in  terms  of  values 


I 


- 
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on  the  preceding  row,  this  formula  is  called  explicit.  If  (4.4)  is 
applied  along  the  row  . . . ,3,0, 1, . . .  we  obtain  a  set  of  simultaneous 
algebraic  equations  which  can  be  solved  to  give  the  solution  along  the 
row  ... ,6,2,5, .. .  .  Equation  (4.4)  is  called  implicit  because, 
although  the  unknown  values  cannot  be  obtained  directly  from  the  equation, 
it  does  imply  a  set  of  equations  from  which  the  solution  may  be  obtained. 

It  is  possible  to  generalize  (4,3)  and  (4.4)  by  taking  their 
weighted  average  as  follows:  Let  6  be  a  constant  such  that 

o  <  e  <  l  , 

Multiply  (4.4)  by  9  ,  (4.3)  by  (l=@)  and  add  to  get 


§4.2,2.  Stability  in  partial  difference  systems: 

Let  us  digress  for  a  moment  and  consider  the  problem  of 
stability  in  partial  difference  equations.  In  §2.4.3  we  defined  in¬ 
stability  in  an  ordinary  difference  equation  as  the  growth  of  parasitic 
solutions  which  eventually  dominated  the  true  solution,  and  we  found 
that  these  parasitic  solutions  were  introduced  by  an  initial  inaccuracy 
in  the  data.  The  properties  of  instability  in  partial  and  ordinary 
systems  are  quite  similar.  Boundary  value  problems  (ordinary  and 
partial)  are  generally  stable,  while  initial  value  problems  are  subject 
to  instability.  Instability  in  partial  systems  generally  manifests 
itself  by  a  rapid  oscillation  of  the  numerical  solution  as  well  as  an 
iircreasing  departure  from  the  analytic  solution  as  one  proceeds  forward 
in  the  time  direction.  Reduction  of  the  interval  in  the  Redirection 
(h  )  without  a  corresponding  reduction  of  the  time  interval  (h  ) 
only  magnifies  the  error.  It  has  been  observed  that  a  drastic 
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reduction  of  the  time  interval  reduces  the  effects  of  instability  and 

2 

that  the  system  is  stable  for  some  values  of  h  /  h  and  instable  for 

x  '  t 

others . 


Quantitative  investigations  of  instability  in  simple  initial 
systems  have  been  carried  out  [1915»16,1Y3.  Hildebrand  [15]  defines 
instability  as  the  tendency  for  the  effect  of  an  initial  numerical 
inaccuracy  (such  as  round-off  error)  to  increase  unboundedly  as  the 
solution  is  carried  forward  in  the  t-directlon. 

§b„2.3°  Stability  in  parabolic  systems: 

Among  the  most  complete  of  stability  analyses  and  of  particular 
interest  to  us  are  the  results  of  Courant,  Friedrichs  and  Lewy  [16] 
relating  to  the  parabolic  equation  (b.2)  with  data 

U(x,o)  =  qp(x)  , 

U(o,t)  S3  U(jc,t)  =  0 

2 

and  a  a  constant.  It  was  found  that  the  difference  equation  (4.6) 
is  stable  if 

2  lit  1  1 

2  a  7 2  «  rfe  •  0  s « <  5  • 

x 

and  with  no  restriction  if 


\  <  B  <  1  . 


The  argument  used  to  obtain  this  result  is  tedious  and  will  not  be 
given  here. 


■  3 
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Richtmyer  [14]  argues  that  the  same  conditions  hold  for  the 

more  general  parabolic  equation 


.2 


du  _  _  „ 

5T  *  a  “  +  b 


Bu 


+  c  U 


where  a  ,  b  and  c  are  constants,  if  the  former  is  weakened  to 


0  2 

2  "  h  2  ^  1-20 

x 


0  <  0  <  \ 


By  writing  (4.1)  in  the  form 


dU  K  „4  ^  „  .,3  oU 

™  =  5  U  _  +  20  ^ 

2  4 

one  can  draw  a  comparison  between  a  and  5  ^  and  argue  that  the 
difference  equation 


(M)  (u2-uo)  -  ht[e6^(u5)+(i-@)B^(ij 


is  stable  if 

,  h 

10  u4  —2  <  .  0  <  e  <  4  , 

X 

with  no  restriction  for  ~  <  0  <  1  .  That  this  is  indeed  the  case 
has  been  verified  by  numerical  work  [14]. 


From  the  form  of  (4.7)  it  is  clear  that,  for  other  than 
0=0,  we  must  solve  non-linear  algebraic  equations.  However,  if 

0=0  the  stability  condition  is 


and. 


unless  U  <  <  1, 


h 

t 


must  be  extremely  small  if  we  are  to  avoid 


a  very  coarse  x  interval. 


. 


Hence,  in  the  numerical  work  below  we  will  employ  the  implicit  scheme 
(0=l)  and  linearize  the  difference  equation. 


Analytic  Properties  of  the  Equation  •  =  «$■“•  ; 

2  at 

OX 

§U„3°1»  Analytic  solution: 

Writing 


U(x, t)  =  X(x)  T(t) 


and  substituting  in  (4.1)  we  obtain 


^  ~r  -  <« 

dx 


Hence, 


1  d2X^ 
Xdx2 


1  dr 

J?  dt 


where  A  is  a  constant, 


Consider  first  the  equation  for  T 


L  dT 
3  dt 


=  -x 


This  integrates  to 


k  T  a  -A  t  -  C 


where  c  is  a  constant. 


Thus 


* 
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If  we  associate  U  with  temperatures,  we  may  assume  that  X 
and  T  are  positive.  We  will  require  both  X  and  T  to  be  real  and 
bounded.  This  means  that  both  A  and  c  must  be  positive.  Thus  the 
t-dependent  part  of  U  is 


(b.8)  T  =  [l*(c  +  At)] 


-b 


where  c  and  A  are  positive  constants, 


The  equation  for  X  is 
“AX, 


i£ 

2 


dx 


where  X  is  positive.  Let 


Y  =  X^  >  0  . 


Then 


dfy 

dx2 


A  Y 


1/5 


Multiplying  both  sides  of  the  above  by  and  noting  that 

[l  (Sfl 


eft  dY_  d_ 

,2  dx  _  dx 
dx 


we  obtain 


d_ 

dx 


[i  @7] 


A  yl/5  « 

X  dx 


After  integrating  with  respect  to  x,  we  have 


(sh 


6/5 


const.  -  ,  A  Y 


fMa-  Y6/5)  . 

5 


where  a  is  a  constant.  Since 


T  >  0  ,  A  >  0  , 


then 


a  >  0  . 


Since  Y  is  bounded  then 


Hence, 


lim  Y 
x— ►go 


6/5 


Let 

Y  =  S5/2  . 


Then 

dY  _  £  s3/2  dS 
dx  2  dx  ’ 


Hence,  the  last  equation  in  Y  can  be  written 


or 


|s3/2i  -  [|X(..s3)jk, 

AA'i  _  SV2  dS 

VV  =  (a-S3)i  dx 


s2  dS 

(aS-S^)^  dx 


If  we  integrate  with  respect  to  x  we  obtain 
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where  b 


and 


are  constants. 


By  letting 


a  =  Qt\3 


we  can  write  this  last  equation  in  the  following  form: 

u  2 


(^•9) 


where 


(0 


x  +  b  =  T] 


r  udu 

J  <8u-i?y^  ’ 


u  = 


[X(x)j‘ 


It  can  be  shown  that  the  integral 

x  =  P  u2du 

V  (8u-UU)2 

reduces  to  a  sum  of  elliptic  integrals  and  the  details  of  this  reduction 
appear  in  Appendix  I.  It  was  originally  intended  that  values  of  the 
analytic  solution  be  compared  with  the  numerical  solution  in  order  to 
illustrate  the  methods  developed  in  this  thesis.  This  would  involve 
numerical  evaluation  of  the  integral  in  (4.9)  or  its  equivalent  elliptic 
integral  sum.  From  a  practical  point  of  view  the  former  is  preferable 
since  tables  of  elliptic  integrals  of  the  third  kind  are  incomplete  and 
interpolation  in  tables  (where  available)  would  be  necessary  in  any 
case. 


§4.3.2.  Discussion  of  the  analytic  solution: 

For  the  purposes  of  this  section  let  us  write  equations  (4.8) 
and  (4.9)  in  a  slightly  different  form  as  follows: 


T(t)  = 


(d  +M,:)'lA 


) 


x  +  b 


n 


X2/>| 

2. 

u  du 

. . ri 

(8u-U  )2 


9 


(^.10) 


- 
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U(x,t)  =  X(x)  T(t) 
is  the  solution  to  the  P.D.E.  (k.l). 

As  pointed  out  in  Chapter  III,  a  P.D.  problem  is  meaningful 
(in  a  practical  sense)  only  if  auxiliary  conditions  are  imposed  on  the 
solution.  Usually  we  must  turn  to  the  physics  of  the  problem  for 
appropriate  conditions  for  it  is  not  wise  to  apply  arbitrary  data  since 
this  often  results  in  a  system  which  is  not  well-posed.  A  physical 
analogue  to  (l+.l)  has  not  been  discovered.  In  any  case,  the  implicit 
nature  of  the  second  of  (k.10)  makes  the  application  of  auxiliary 
conditions  difficult.  Therefore,  let  us  reverse  the  normal  procedure 
and  examine  the  solution  so  as  to  discover  what  data  is  appropriate 
to  the  equation. 


We  have  at  our  disposal  the  constants  d,  b,  q ,  u^  and  |i. 
The  origin  in  x  and  t  is  determined  by  d  and  b,  for 


T(o)  =  d"1/^ 


X2(o)/q 


2, 

u  du 


b 


u 


1 


It  is  reasonable  to  take 


T(o)  =  1 


that  is 


d  =  1  . 
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Since  is  as  yet  undetermined,  we  can  let 

b  =  0  . 

The  rate  at  which  the  solution  decays  as  t  increases  is 
determined  by  4.  It  is  convenient  to  set  |i  =  15.  This  makes  T  decay 
fairly  rapidly,  but  a  full  investigation  would  require  a  range  of 
values . 


The  initial  shape  of  the  solution,  U(x,o)  ,  is,  in  the  main, 
fixed  by  T|.  Again,  a  number  of  values  should  be  tried,  but  tj  =  1 
is  a  reasonable  value  to  investigate. 


Our  solution  now  has  the  form 


T  =  (1  +  150"1/4  , 


(4.10a) 


x  = 


2, 

u  du 

FT 

(8u-u4)2 


We  have  required  that  X  be  a  real  function  of  x,  hence, 
we  must  have 

k 

8u  -  u  >  0  . 

Now 

k  p 

8u-u  =  u(2-u)(u  +2u+4)  , 

but 

u2+2u+4  =  [(u+l)2+3] 

which  is  always  positive.  Hence,  we  must  have 


0  <  u  <  2  . 
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Let 


It  will  now  be  shown  that  X  is  a  periodic  function  of  x, 


w  =  u  +  iv 


be  a  complex  variable  and  consider  the  integral 

2 


P  w  d 
J  ( 8w-w 


dw 

k  JL.  » 

(8w-w  )2 


where  C  is  a  "dumbbell"  shaped  contour  composed  of  the  circles, 
and  C^,  about  the  branch  points  w  =  2  and  w  =  0  joined  by 
the  lines  AB  and  CD  as  shown  in  the  figure  below. 


Suppose  that  the  circles  have  radii  e  and  that  the  contour  is  so 
chosen  that  all  other  singularities  are  excluded.  We  may  write 


r  w^dw  r  w^dw  P  w^dw  r  w^dw 

^  ( 8w“W^) 2  (8w“W^)2  ^  (8w-w^)2  (8W-W11)2 


2 

1W 

(8w“W^)‘ 


2 

J  c 

(Sw-W^)' 


On  we  write 


2  -  w  =  e  e 


icp 


Then 


r .  - 1  «*  r  (2-^1 

\J  J  , ,  _  i®v 


iq>.3/2oi<p/2 


dcp 


it  [(2-cei<p)S+2(2-ee1<P)+^ji 


and  it  is  clear  that  this  integral  vanishes  when  €  -+■  0  .  Note  that 


the  function 
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2 

 w 

[w(w  +2w+l+)]2 

is  unchanged  in  sign  as  w  passes  around  „ 

Now 

(2-w)^  =  e2  eiq>/2  . 

When  w  is  at  A,  cp  =  -  jt  and 
(2-w)2  =  -  i  . 

When  w  is  at  C,  cp  =  jt  and 

(2-w)2  =  i  . 
l 

Hence,  (2-w)2  changes  sign  as  w  goes  around 
similar  arguments  for  w  on  the  circle  C^. 
Thus,  if  we  let 


v  -*■  0,  e  0  , 


w  dw 
—  TTi 
(8w-w4)2 


2, 

u  du 

- -  > 

(8u-u4)2 


and 


f  w2dw 
CD  (8w-w^)^ 


u  du 
(8u-u^)^ 


I  = 


u  du 
(8u-u^ 


Hence,  as  u  goes  through  the  range  0  to  2  to  0, 


2, 

u  du 
'  h~±. 

(8u-U4)2 


.  We  can  apply 


x  increases  by 


2 
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Therefore,  u  is  a  periodic  function  of  x  with  period 
r 2  u2du 

2  /  - - . 

(8u-u4)2 

S  ince 

X(x)  =  [u(x)]2  , 

then  X  is  a  periodic  function  of  x  . 

The  point  of  the  preceding  analysis  is  this:  just  as  we 
can  investigate  the  function 


y  =  sin  x 


by  considering  the  integral 


ry  at 

*  =  l  d-t2)4 


for 


l  <  y  <  l  » 


0  <  x  <  2jt  =  2  I  - ~i  , 


r1 _ at, 

J,  (i-t2: 


we  can  discover  the  properties  of  X(x)  by  considering 

,„2 


x 


for 


2, 

u  du 
(8u-uU)i 


0  <  X  <  2 


r2  2, 

/  u  du 

Jo  (8U-U21)2 


0  <  x  <  2 
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We  can  complete  the  problem  by  evaluating  the  above  integrals 
numerically. 


A  Numerical  Procedure  for  Solving  the  Equation 

§  1+  - 4 . 1  -  Initial  conditions; 

In  view  of  a  reasonable  data  curve  to 


d2U^  dU 


choose  would  be 


t  =  0, 


9 


where 


x 


2 


> 


The  initial  conditions  are 


U(x,o)  =  A(x)  , 
U(xrt)  =  B(t)  , 
U(x2,t)  =  C(t) 

and  the  region  of  solution  is 


t  >  0  ,  <  x  <  x^  . 


§ 4 o 4 . 2.  The  difference  equation: 

To  avoid  difficulties  arising  from  instability  let  us  use  an 
implicit  scheme  ( 8  =  1  in  equation  (^.7))-  The  difference  equation  is 


k  52  u5.  -  Vr 

X  1J  t 


U(x1+ihx,jht)  =  U. .  , 

k  =  V^X  > 


where 


-  6o  - 


il  ■  k  <5x-6x+---)  Uii  +  <fVt+^Vt  +  • 


-j 


ij 


i  ”  1 1  2  j  •  •  «  ,  n  ^ 
J  ~  I?  2j  » e « j  nij 


x  -x 

h  =  h  =  ~ 

x  n+ 1  t  m+  1 


and 


U,  =  A  ,  U  .  =  B  . ,  U  _.  =  G. 
io  i  oj  j  n+l,j  j 


We  must  now  linearize  the  difference  equation.  This  can  be 
done  in  several  ways.  First,  let  us  write  the  equation  explicitly  in 
terms  of  the  unknown  quantities, 


u.  .  u, U,  .  . 

1-lfJ  ij  i+l»J 


as  follows: 


k  .  ,  -  2k  IT?.  -  U. .  +  k  uj  .  .  =  -U.  .  -  +  c  . 
i“l»J  ij  ij  i+l,J  i»J-l  ij 

4 

If  we  let  be  an  approximation  to  U  in  terms  of  known  quantities, 

we  have 


ka.  ,  U.  ,  .  +  (3.  U.  .  +  ka.  ,  U.  .  .  =  -U.  .  -  +  c,  . 
i-1  i“l,j  'i  ij  i+l  i+l»  J  isJ“l  ij 


where 


|3.  =  -(2k  +  l) 


The  difference  equation  in  matrix  form  is 


(4.11)  FZ  =  G  , 


where 
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Pi 

ka2 

0 

0 

o 

0  0 

ka^ 

P2 

ka, 

j 

0 

9 

0  0 

F  = 

0 

ka2 

ka, 

4 

9 

0  0 

0 

0 

e  o 

0 

0 

« 

9  9  9 

.  ka  ,  P 

n-1  _n. 

Ulj 

C1  . 
Ij 

-  U-  .  ,  -  kBl 
l.J-1  J 

Z  = 

U2j 

c 

3 

G  = 

'<—1 

OJ 

o 

-  U„  .  - 

2,  J“1 

o 

u  . 

c  . 

__nJ 

4 

-  U  .  ,  -  k€7 

n,J“l  J_ 

a 

o 

and 

an+l 

are  known  exactly, 

that  is. 

a  = 

o 

4 

U  . 
OJ 

=  B 

J 

3 

4  4 

a  .  =  IT  -  .  =  C7 
n+1  n+lsj  j 


One  linear izationj  based  on  the  assumption  that  U  is  a  slowly  varying 
function  of  y,  is  to  let 


I. 


ai  ~  Ui,j“l 


We  then  compute  the  vector  Z  thereby  obtaining  a  first  approximation, 

( l) 

'  ■'U  .  .  Next  we  let 
ij 

„  _  (i)„ 

ai  “  Uij 


and  solve  again  to  get  a  second  and  (we  hope.)  improved  approximation, 

^2V..  With 
13 

„  _  (2)„ 

ai  "  “ij 
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we  solve  (4.11)  again,  repeating  the  procedure  until  successive 
solutions  agree  to  the  required  accuracy.  Initially  we  must  let 
Cij  =  but  dif ference-corrections  can  be  applied  at  any  stage  after 
the  first. 

Since  the  explicit  difference  equation  is  linear  we  can 

..  If  we  then  let 
ij 


easily  calculate 


„  _  (Du 

ai  “  uij  ’ 


we  may  be  able  to  improve  this  first  estimate  by  iterating  as  before. 
Although  the  explicit  equation  would  probably  be  instable,  the  error 
should  not  be  too  large  for  only  one  time  step. 

Richtmyer  [14]  linearizes  as  follows:  using  Taylor's  theorem, 


U?.  =  B5 
ij 


i,j-l  +  ht  )u.^  4 


/ 


T  «  «  o 


,t5  cit4  .  f  du 

=  U,  .  ,  +  5U.  .  ,  h  ( 

i, j“l  i , j=l  t  y  ot 


i.J-1 


but 


h  (P  )  „  A  »,  .  , 

~  c  i*j-x 


Hence, 


U.  .  -  U.  .  . 
~  i,j  i,j“l 


u?.  u?  .  .  +  5U1!  .  ,  (u.  .-u.  .  ' 
ij  ~  i,j-i  i>j~i  i»j  i»j-i 


Vj  -  • 


-  63  - 


§4,4.3*  Some  numerical  results; 

An  algorithm  was  written  to  solve  the  following  system: 


S2u5  Su 


(U.12)  u(0,t)  =  0(1, t)  =  0  , 

U(x,o)  =  sin  jt  x  , 

0  <  x  <  1  ,  t  >  0  . 

The  differential  system  was  reduced  to  the  difference  system  (4.11) 

and  the  first  linearization  scheme  of  §4.4.2  was  used.  The  first 

(l) 

approximation  to  the  first  row  solution,  '  'U  ^  s  and  three  success¬ 
ive  approximations,  are  given  in  Table  II. 

The  mesh  size  is 

hx  -  ht  -  Vto  • 

The  solution  is  symmetric  about  x  =  0.5  and  only  half  of  the 
solution  is  recorded. 

TABUE  II 


Approximate  solution  of  the  system  (4.12)  along  the 
x  ^U(x,l/10)  ^U(x,l/10)  ^0(x,l/lO) 

line  t  =  1/3 
^U(x9l/10) 

0 

0 . 0000 

0 . 0000 

0 . 0000 

0.0000 

0,1 

1.1271 

0.0317 

0.3165 

0.6881 

0.2 

0.8576 

0.1377 

2.0933 

0.0269 

0.3 

0.5169 

0.7398 

0 

0 

— >3 
O'! 

0.5729 

0.4 

0.3805 

1.1569 

0.1525 

1.3261 

o.  5 

0.3441 

1.1598 

0.1738 

0.9961 

' 
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The  successive  solutions  show  no  tendency  to  converge. 

Decreasing  the  mesh  size  effects  some  change,  but  mesh  sizes  of  the 
order  of  1/100  would  be  necessary  for  significant  improvement. 

Clearly  the  linearization  procedure  used  is  too  crude.  It  is  likely 
that  the  choice  of  initial  data  is  inappropriate,  introducing  a 
rapidly  varying  dependence  on  t  into  the  solution. 

Consideration  of  this  problem  is  far  from  complete,  but  the 
foregoing  has  served  to  point  out  the  inadvisability  of  considering 
the  differential  equation  as  a  mathematical  abstraction  with  arbitrary 
initial  conditions.  Lacking  a  physical  basis  for  &  meaningful  problem 
suitable  data  will  have  to  be  obtained  numerically  from  the  analytic 
solution  using  various  values  for  the  constants 5  only  then  will  it  be 
possible  to  undertake  the  numerical  solution. 
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CHAPTER  V 

ELLIPTIC  SYSTEMS 


§5ol.  Introduction: 

The  wide  application  of  closed,  difference  methods  to  the 
solution  of  boundary  value  problems  in  ordinary  differential  equations 
[5]  would  lead  one  to  suspect  that  similar  methods  are  commonly  used 
in  the  solution  of  elliptic  systems.  This  is  not  the  case.  To  date, 
matrix  methods  have  been  applied  only  to  a  limited  class  of  linear 
equations  with  very  simple  boundaries,  whereas  relaxation  procedures 
[3]  and  sweeping  techniques  [20]  Recent  improvements  of  relaxation) 
have  received  more  attention  and  have  been  successfully  adapted  to  many 
types  of  equations  and  boundaries.  Nevertheless,  matrix  methods  offer 
a  conciseness  of  formulation  and  an  invariance  of  procedure  which  are 
lacking  in  sweeping  techniques. 

In  this  chapter  some  of  the  matrix-difference  methods  for 
solving  linear  elliptic  P.D.  systems  are  discussed.  At  first  only  very 
simple  equations  are  considered,  but  generalizations  are  made  wherever 
possible.  We  assume  that  the  boundary  curve  is  a  rectangle  with  sides 
parallel  to  the  axes  and  that  the  solution  is  specified  on  the  boundary. 
Non-rectangular  boundaries  and  boundary  conditions  involving  the  normal 
derivative  of  the  solution  are  considered  separately.  One  section  of  this 
chapter  is  devoted  to  a  discussion  of  singularities  in  boundary  conditions, 
but  otherwise  it  is  assumed  that  the  solution  and  data  are  continuous  and 
non-singular.  In  most  cases,  concepts  are  illustrated  by  particular 

1 

examples.  In  what  follows  we  will  use  results  from  the  theory  of  matrices 
[6,18];  no  proofs  for  these  theorems  will  be  given. 


-  66  - 


§5.1.1.  Conventions; 

To  avoid  unnecessary  repetition  in  later  work  let  us  establish 
the  following  conventions: 


The  region  of  solution  of  the  equation  is  the  rectangle 


a  <  x  <  a.  ,  b  <  y  <  b, 
o  “l  o  1 


A  grid  or  mesh  is  set  up  in  this  region  consisting  of  the  points 


(*.  ,  y^)  , 


where 


x.  =  a  +  ih  , 
l  o  x 


y .  =  b  +  jh  , 

J  o  y 


i  =  0,1,... ,n+l. 


j  =  0,1,... ,m+l , 


a  -a  b  -b 

.  _  _1 _ o  .  _  _1 _ o 

x  ”  n+1  '  y  m+1 


Interior  points  are  those  points,  (x^,y^),  for  which 


x  —  l,2,o..,n, 

j  —  1,2, ... , tn 


and  boundary  points  are 

(xi’yo)>  (xi'Vl>  1 

X  =  0  j  1  )  »  ♦  «  )  H-j*  1  j 
j  =  0,1,... ,m+l  . 


- 
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The  value  of  the  solution  at  the  point  (x^,y^)  is  denoted 
by  z-jj*  The  matrix  of  the  solution  at  interior  points  is  denoted  by 
Z  or  (z_).  We  also  use  the  letter  Z  to  represent  the  contin¬ 

uous  solution  whenever  a  P.D.  system  is  written  down.  It  should  be 
clear  from  the  context  whether  Z  stands  for  a  matrix  of  discrete 
values  or  the  continuous  solution. 


Other  conventions  will  be  given  as  we  need  them. 


§5.2.  Matrix  Representations  of  Some  Simple  Elliptic  Systems: 

Among  the  most  common  of  elliptic  systems  are  those  involving 
Poisson's  equation.  The  method  used  to  reduce  this  equation  to  a  matrix 
system  can  easily  be  extended  to  handle  more  general  types  of  elliptic 
P.D.E's  and,  hence,  will  serve  as  a  model  for  later  work. 


(5.1) 


Let  us  write  Poisson's  equation  in  the  form 

^-§  )  =  f(x.y). 

Sv2  / 


sfz  a^z 
a*2  +  ay‘ 


The  reason  for  the  minus  sign  will  be  apparent  later.  Multiply  (5.1) 

2  2 

by  (h^  hy)  and  use  the  central  difference  formula  (2.16)  to  obtain 
the  difference  equation 

(5.2)  -  (h2S2z, .+h262z, .)  =  h2h2f  +(h2c  +h2c  )z, .  , 

w  '  y  x  ij  x  y  i]'  x  y  ij  '  y  x  x  y'  ij 

where  c  and  c  are  difference-correction  terms, 
x  y 


It  is  established  custom  to  have 


h  =  h  =  h 
x  y 


unless  it  is  known  in  advance  or  by  a  trial  computation  that  the 
solution  varies  more  rapidly  in  one  direction  than  in  the  other. 
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Let 


(c  +  c  )z .  .  =  c .  .  , 
x  y  ij  ij 


h  fij  +  Cij  “  8ij  * 


Then  (5-2)  can  be  written 


(5-3) 


(6xZij+8yZlj>  -  *ij 


i  —  1)2) >  * « ^n  , 
j  =  1 , 2, . . . ,m  . 


A  double  array  of  quantities  such  as  (5. 3)  can  be  most  conveniently 
represented  by  a  matrix.  Let  us  consider  the  special  case  n  =  3, 


m  =  4  to  see  if  a  compact  matrix 
found.  In  matrix  form,  therefore 


-zoi+sairz2i 

‘Z02+2z12"Z22 

■Zll+2Z2l'Z31 

-Z12+2Z22-Z32 

‘Z21+2z31'ZU 

-Z22+2Z32-V 

~zio+2zirzi2 

-ZU+2Z12-Z13 

"Z20+2Z21"Z22 

"Z21+2Z22'Z23 

^z30+2z31'z32 

"z31+2z32"Z33 

11 

S12 

813 

8U 

21 

§22 

g23 

824 

_31 

g32 

S33 

83U 

Now,  the  two  matrices  on  the  left 
contain  the  given  boundary  data, 


representation  of  (5*3)  can  t>e 
we  have 


~z03+2z13~z23 

-Z04+2z1U-Z24 

-213+2Z23'Z33 

-zih+2zsk-z}k 

-z23+2z33"Z43 

-*2h+2z}h-Zkk 

-Z12+2Z13'Z14 

-z13+2z14'z15 

"Z22+2z23"Z2lv 

‘Z23+2z24"Z25 

•Z32+2Z33*Z34 

•Z33+2z34-z35 

hand  side  of  the  equation  above 
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Z  .  •  Zl  ,  •  Z  )  Z  y 

Oj'  4j  io’  i5 


i  =  1,2, 3, 4, 
j  =  1,2, 3, 4, 5. 


Let  us  retain  only  unknown  quantities  on  the  left  side  of  the  equation 
by  adding  the  matrix 


or  io  z02  03  o4+  15 

0  0 


20 


25 


L_!41+z30  z42  z43  Z44+z^5 


to  both  sides c  The  matrix  equation  is  now 


z2  =  F  , 


where 


2zll"Z21 

2Z12"Z22 

2z13“Z23 

2z14"z24 

X  = 

-zn+2zsi’z3i 

-Z12+2Z22-Z32 

"213+£223'Z33 

-Zl!*+2z2!t- 

2  3^ 

'z21+ 2z31 

"Z22+2Z32 

•Z23+2Z33 

■Z2l*+2Z31* 

2z11"Z12 

"ZU+2Z12-Z13 

"Z12+2Z13"Z14 

•Z13+2ZuT 

II 

CVI 

fr* 

2z21"Z22 

•Z21+2Z22'Z23 

-Z22+2Z23'Z21. 

-Z23+2Z2lt 

2z31"z32 

'z31+2z32'Z33 

•Z32+2Z33‘Z3^ 

•z33+2z3_^ 

gn+zoi+zio 

g12+Z02 

g13+z03 

814+z04+z15 

F  = 

g2l+Z20 

g22 

823 

g24+Z25 

831+z41+z30 

832+z42 

833+Z43 

g54+z44+z35 

The  matrix  can  be  written  as  the  product 


:o  &b:  &  :  aril  no  sJ-  3- •  -p  nwot  v  ^  V  zv  ?*j. 
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~2 

-1 

o- 

Z11 

Z12 

213 

zlb 

-1 

2 

-1 

Z21 

Z22 

Z23 

Z2b 

0 

-1 

2 

_ 

Z32 

Z33 

N 

1  ^ 

and 


y  can  be  written  as  the  product 


11 

z12  Z13 

zlb 

21 

Z22  Z25 

z2b 

31 

Z32  z33 

Z3b 

2-100 
-1  2-1  0 

0-12-1 
0  0-1  2 


Following  Bickley  and  McNamee  [19],  let 


where  the  subscript  denotes  the  order  of  the  matrix.  Thus,  we  can 
write  the  matrix  equation  as 

2  2 

Df  Z  +  Z  Df  =  F  . 

3  ^ 

2 

Let  be  the  triple  diagonal  matrix  of  order  (nxn)  whose 

diagonal  elements  are  2  and  whose  sub-and  super-diagonal  elements  are 
-1.  Then,  by  the  same  procedure  as  employed  in  the  example  above,  we 
can  write  (5-3)  in  matrix  form 

(5.M  D2  Z  +  Z  D2  =  F  , 

'  7  n  m 
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where  Z  and  F  are  matrices  of  order  (nxm).  F  consists  of  the 
elements 

(h2fi.  +  Cij}  , 

with  the  quantities 

oj  n+l,j  io  i,m+l 

added  to  the  first  row,  last  row,  first  column  and  last  column  respect¬ 
ively,  where 

x  —  l,2,ooo, n , 

J  “  l,2,ooo, m o 


It  should  be  pointed  out  that  equation  (5.U)  represents  the 

complete  difference  system  (equation  with  boundary  values).  We  also 

2 

note  that  D  is  a  positive-definite  matrix.  It  was  to  assure  the 
n 

positive-definiteness  of  that  we  chose  the  minus  sign  in  (5.1). 

In  precisely  the  same  way  we  can  represent  the  equations 


(5-5)  "  (c(x)  “-|  +  e(y)  7“|  ^  =  f(x,y) 


6x£ 


and 


(5-6) 

by 


d2[e(x)z]  6fXeXz)Zl 


dx‘ 


^  =  f(x,y) 


( 5 . Ua)  CD2  Z  +  Z  D2  E  =  F 
v  '  n  m 


and 


(5,^b)  D2  CZ  +  ZED2  =  F 
'  '  n  m 
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respectively.,  where  C  and  E  are  diagonal  matrices  of  order  (nxn) 

and  (mxm)  respectively „ 

§5,3,  Solution  of  the  Matrix-Difference  Equation; 

Equations  (5„lf),  (5Aa),  (5„4b)  can  be  written  in  the  form 

(5.7)  AZ  +  ZB  .  F  , 

where  A  is  of  order  (nxn),  B  is  of  order  (mxm)  and  Z  and  F 

are  of  order  (nxm),  A,  B  and  F  are  known. 

Three  methods  of  obtaining  the  solution  of  (5. 7)  were  given 
by  Bickley  and  McNamee  [19].  The  first  of  these  involves  the  eigen- 
quantities*  of  both  A  and  B  and  is  called  the  irrational  method. 

The  second,  termed  semi-rational,  requires  the  eigenq  > title--,  of  A 
or  B  but  not  both.  The  third,  named  rational,  does  not,  in  theory, 
require  eigen quantities  at  alL. 

If  a  matrix  of  order  (kxk)  has  k  linearly  independent 
eigenvectors 9  it  is  defined  to  be  non-defective ,  In  the  irrational 
solution  both  A  and  B  must  be  non-defective,  while  the  semi- 
rational  method  requires  that  at  least  one  of  them  be  non-defective. 

In  the  third  method  no  restriction  is  placed  on  A  or  B.  The  reason 
for  the  names  irrational,  semi-rational  and  rational  will  become  appar¬ 
ent  later, 

§5,3,1,  The  irrational  solution; 

Suppose  that  A  and  B  are  non-defective  and  that  A  has 
the  n  eigenvalues  X  ,  the  n  eigenrows  and  the  n  eigen- 


* 


Eigenquantities-eigenvalues  and  eigenvectors. 


-  73  - 


columns  I\  while  B  has  the  m  eigenvalues  v ^ ,  the  m  eigen- 

rows  cr.  and  the  m  eigencolumns  s.. 

J  J 

The  eigenquantities  of  A  obey  the  following  relations: 


°t  A  =  \  pt  > 


(5.8)  A  r.  =  X.  rt  , 


p .  r .  =  5. . 

1  J  ij 


The  eigenquantities  of  B  obey  the  relations 


cr.B  =  v.  cr. 
J  J  J 


9 


(5.9)  Bs.  =  v.  s.  , 


cr.  s . 
1  1 


5 


ij  ° 


Let  u. .  be  scalars,  where 
ij 


i  =  1,2, 

j  =  1  S  2  s 


*n  , 

5  m  o 


The  theory  of  matrices  tells  us  that  an  arbitrary  matrix  can  be  repres¬ 
ented  as  a  linear  combination  of  products  of  the  eigenvectors  of  non¬ 
defective  matrices „  In  particular 

n  m 


Z  = 


I  I 

i-1  i-1 


r.  cr.  . 
1  J 


If  we  substitute  this  in  (5.7)  we  have 


A 


I  uij  ri  + 

i»  j 


i»  j 


riV  =  F‘ 


In  view  of  (5.8)  and  (5.9)  this  becomes 


I  (V  uij  ri  ■ F  • 


Now,  if  we  pre-multiply  this  relation  by  p^,  post-multiply  by  s 
and  use  the  last  of  (5.8)  and  (5*9) »  we  obtain 


(\  +  VP  "u  =  °k F  si  • 

Hence,  provided 


(\  +  vj)  /  0  > 


we  can  find  the  scalars  u, . 

ij 


p .  Fs  . 

(A7+vTi  » 


from 


i  —  1,2, ... ,n  , 

j  "■  l,2,ooo, tn  . 

This  method  requires  that  we  find  the  eigenquantities  of 

both  A  and  B,  which  is  no  mean  task.  So,  except  for  a  few  special 

2 

cases  (in  particular  D^),  this  method  is  to  be  avoided  in  practical 
computation. 


§5.3.2.  The  semi°rational  solution: 

Let 

Uj  ,  j  =  1,2, ... ,m 

be  columns  of  order  (nxl).  Suppose  that  B  is  non-defective,  but 
place  no  such  restriction  on  A.  We  can  expand  the  solution  matrix, 
Z,  in  the  form 
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CTj  B  =  F. 

(5.9)  we  have 

U .  <7 .  =  F  . 

J  J 

last  of  (5*9)  to  get 

Therefore,  if  we  define  I  to  be  the  matrix  with  unit  diagonal  and  zero 
off-diagonal,  U.  can  be  obtained  by  solving  the  equation 

(A  +  v.l)  U  =  F  s.  , 

1  i  1 

provided 

|A  +  v.l|  /  0  . 

This  is  equivalent  to  the  condition 
+  v±)  /  0  . 

It  is  clear  that  a  similar  relation  can  be  found  in  terms 
of  the  eigenquantities  of  the  (non-defective)  matrix  A. 

The  semi-rational  solution  is  of  more  practical  value  than 
the  irrational  solution  even  though  n  sets  of  m  linear  equations 
in  m  unknowns  must  be  solved,  for  obtaining  eigenquantities  accur¬ 
ately  is  exceedingly  tedious. 


Z  =  )  U.  <7. 

L j  j 
j=i 


Equation  (5.7)  becomes 


A  )  U.  or.  +  )  U. 

Z_j  j  j  Z_i  j 


After  application  of  the  first  of 


A  )  U.  a.  +  )  v. 

Z_j  j  j  Z.J  j 


Post-multiply  by  s^  and  use  the 


A  U„  +  v.  U.  =  F  s  . 
111  i 


33§  GJ  (^.?)  5 o  jbJ  sd:>  9;  '  :U> 
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§5.3*3*  The  rational  solution: 

Before  proceeding  with  the  rational  solution  of  (5*7) »  a  few 
results  from  matrix  theory  should  be  stated.  We  defined  the  eigen- 
quantities  of  a  matrix  B  such  that  they  satisfy  the  relations 


a,  B  =  v.  o\  , 
i  l  i 


B  si  -  vt  si  • 


These  equations  can  be  written  as 


°’i(B  -  v.l)  =  0  , 
(B  -  Vil)  s£  =  0 

and  hence,  the  determinant 


| B  -  v^l  =  0  . 

Expanding  the  determinant  gives  us  a  polynomial  in  v  which  is  called 
the  characteristic  polynomial  of  the  matrix  B  and  is  written 


m  m-1 

v  -  p^v  + 

(5.10)' 

or  (v-vx)(v-v2)  . 


+  (-D-p.  -  o 

(v-v  )  =  0. 

'  m 


The  Cayley-Hamilton  theorem  states  that  a  matrix  satisfies  its  own 
characteristic  equation,  that  is, 

m  __m-l  /  i t  ^ 

B  -  p^IB  +  .  .  .4-  (“1)  P^I  =  ^ 

(5.10a) 

or  (B-Vjl)(B-v2l)  ...  (B-vml)  =  0  , 


where  0  is  the  zero  matrix. 


■ 
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The  rational  solution  proceeds  as  follows:  let  be 

any  constant.  Add  v^IZ  to  both  sides  of  (5.7)  to  obtain 

(5.11)  (A+Vl)Z  =  F  +  Z  (Vl-B)  , 

where  (A+v^)  implies  (A+v-^l).  Let  be  a  constant  and  multiply 

(5.11)  by  (A+Vg)  to  get 

(5.12)  (A+v2)(A+Vl)Z  =  AF  +  v  F  +  v2Z(Vl-B)  +  V]AZ  -  AZB  . 

But  from  (5.7) 

AZ  =  F  -  ZB  . 


Hence , 


and 


AZB  =  FB  -  ZB 


v^AZ  =  VjF  -  v^ZB  . 


Substituting  these  last  relations  in  (5.12)  we  have 


(a+v2)(A+v1)Z  =  (vj+v2)F  +  (AF-FB)  +  Z[v2(v1-B)  -  V]B+B2] 

=  (v1+v2)F  +  (AF-FB)  +  Z(Vl-B)(v2-B)  . 

Multiplying. this  by  (A  +  v  )  and  proceeding  as  before  we  obtain 

(5.13)  (A+v^)(A+v2)(A+Vl)z  =  (v1v2+v2v^+v^v1  )F  +  ( v^v^v^)  (AF-FB) 

+(a2f-afb+fb2)  +  z(V;l-b)(v2-b)(v3-b)  . 


We  can  continue  in  this  way  with  constants  v,  ,  v..,  ...»  v  ,  but  to 

J  4  5  m 

see  more  clearly  where  this  analysis  is  leading  consider  the  special 
case  m  =  3.  From  (5.10a)  we  have 
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(Vl-B)(v2-B)(v3-B)  =  0  . 

From  (5.10)  we  see  that  the  coefficient  of  F  in  (5.I3)  is  p and 
that  the  coefficient  of  (AF-FB)  is  p^.  Now,  let 

F(k)  =  A1^  -  Ak_1FB  +  Ak~2FB2  -  • • •+  (-l)k  FBk  . 

We  can  write  (5.I3)  in  the  form 

(A+v3)(A+v2)(A+v1)Z  =  p2  F(0)  +  px  F ( 1 )  +  F(2)  . 

For  general  m  ,  it  follows  that 

m-1 

(A+vnl)^A+vm_1)  *  *  '  (A+vl)Z  =  Y  Vk-1  F(k)  ' 


(A+v  )(A+v  -  )  •  •  •  (A+v.  )  =  C  , 

nr  m-17  l7 

the  solution  for  Z  is 

m-1 

z  - C'1  I  Vk-1  F'k)  • 

k=o 


Clearly,  we  could  have  obtained  a  similar  result  in  terms 
of  the  matrix  A  and  in  general  we  will  use  the  matrix  of  lower 
order. 


Note  that  the  matrix  C  can  also  be  given  in  terms  of  the 
coefficients  of  the  characteristic  polynomial  and  so  a  knowledge  of 
the  eigenvalues  of  A  or  B  is  not  necessary.  The  coefficients  of 
the  characteristic  equation  can  (theoretically)  be  calculated  by  purely 
rational  operations,  hence  this  method  is  called  rational.  The  eigen- 
quantities  of  a  matrix  cannot  be  obtained  by  rational  means.  There 
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are,  however,  no  good  numerical  methods  for  finding  the  coefficients 
of  the  characteristic  equation  which  involve  appreciably  less  labor 
than  finding  the  eigenvalues  and  the  distinction  between  rational  and 
irrational  has  little  practical  significance.  The  second  method  is 
called  semi-rational  merely  to  distinguish  it  from  the  first. 


$5.4.  Matrix  Representation  of  More  General  Elliptic  Equations: 

In  §5.2  and  §5.3  we  showed  how  a  restricted  class  of  linear 
elliptic  systems  could  be  represented  in  terms  of  matrix  equations. 
The  purpose  of  this  section  is  to  reduce  equations  of  the  form 

(5.14)  -  (a(x,y)  ^-|  +  b (x,y )  ^-§  ^  =  f(x,y) 

\  dx  dy  S 

to  matrix  equations  and  to  put  forward  one  method  for  solving  them, 
provided  a  and  b  are  slowly  varying  functions  of  x  and  y. 


§5.4.1.  Matrix  equation  for  (5.14): 

The  difference  equation  for  (5.14)  is 


(5.15) 


(aij  8x  zij  +  bij  6y  zij>  ■  gij  • 


where 


g.  .  =  h  f .  .  +  correction  terms, 
ij  iJ 


We  define  the  matrix  dot  product,  C'E  =  L  ,  as  follows:  let 


Then 


0  =  {Cjj}.  E  =  {e^},  L  = 


U 


‘ij  -  cij  eij  • 


that  is,  corresponding  elements  of  C  and  E  are  multiplied  together. 
The  result  (L)  is  a  matrix. 
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If  we  write  the  complete  set  of  difference  equations  (5.15) 
in  matrix  form  as  we  did  in  §5-2,  we  obtain  the  matrix  equation 


(5.16)  A  •  (D2Z)  +  (ZD2)  •  B  =  F  , 
'  '  n  '  '  m' 


where 


A  =  {a  }  ,  B  =  {b.  .) 
*■  ijJ  ij 


and 


F  = 


correction  terms) 


incremented  in  the  first  row,  last  row,  first  column  and  last  column  by 


a.  .  z  . ,  a  .  z  .  . ,  b . ,  z .  and  b  .  z .  n 

lj  oj  nj  n+l,j  ll  10  lm  i,m+l 

respectively. 

§5.4.2.  Solution  of  the  matrix  equation  (5.16): 

Equation  (5.16)  is  to  be  solved  for  the  (nxm)  matrix  Z, 
but  the  nature  of  the  dot  product  does  not  permit  direct  application 
of  the  methods  of  § 5 • 3 •  Let  us  examine  the  dot  product  A'Y  for  the 
special  case  n  =  2,  m  =  3. 


Suppose,  first  of  all,  that  A  has  the  special  form 


Then 


all 

all 

ail 

a 

a 

22 

22 

22 

allyll  a 

lly  12 

ally13 

a22y21  a22y22  a22y23 


an  0 

yll  y12  y13 

a22 

y21  y22  y23 
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that  is, 


-Y  -  *D  Y  ’ 


where  is  a  diagonal  matrix  of  order  (2x2).  Now,  if  A  is  of  the 

more  general  form 


A  = 


we  can  write 


where 


Thus 


Ai  - 


A2  ” 


ail  al2 

a13 

a  ,  a 

a 

21  22 

23 

.  +  A.  . 

f— 1 

ro 

V 

a“ 

Lp 

9 

ro 

a2 

airai 

a12"a 

a2l”a2 

a_  -a 
22 

A1'Y  +  A2- 

Y 

A._  Y  +  A. 

•Y  . 

ID  2 

dot  product 

Y-B 

YBid  +  Y.B 

2  * 

where 
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B 


ID 


B2  = 


1  0 

0 

P2 

0 

9 

0 

P3 

ll"pl 

b12- 

P2 

22_P1 

b22' 

P2 

Hence  in  the  case  n  =  2,  m  =  3  equation  (5.I6)  can  be  written  as 

A1DD2Z  +  ZD^1D  =  F  -  [A2'(D2Z)  +  (ZD3)'b2)  • 

This  suggests  the  general  form  to  which  (5.15)  can  be  reduced,  namely: 


(5.17)  A1DDfZ  ♦  ZD2BlD  .  F  -  [A2.(D2Z)+ (ZD2).B  ]  , 


m  ID 


where  is  a  diagonal  matrix  of  order  (nxn)  with  diagonal  elements 


j  ot*^  * 


,  an  , 


is  diagonal  of  order  (mxm)  with  elements 


h*  P2*  Pm  * 


A  is  of  order  (nxm)  whose  ith  row  is  the  ith  row  of  A  decreased 
2 

by  and  B^  is  (nxm)  whose  jth  column  is  the  jth  column  of  B 

decreased  by  (3  . 

If  a(x,y)  and  b(x,y)  are  slowly  varying  functions,  we 
may  be  able  to  choose  the  and  fL  in  such  a  way  that  the  term 


tV(DnZ>  +  (ZD«>'B21 


can  be  considered  to  be  a  small  perturbation.  In  this  case  we  can 
ignore  the  perturbing  term  and  solve 


by  the  methods  of  §5. 3*  Having  obtained  Z  we  can  evaluate  the  per¬ 
turbing  term,  include  it  in  (5.I7)  and  solve  for  a  better  estimate  of 
Z.  Continuing  in  this  manner  we  may  expect  that  the  procedure  will, 
under  certain  conditions,  converge  to  the  required  solution. 

It  must  be  pointed  out  that  this  method  is  as  yet  untried 
and  the  problem  of  convergence  and  the  selection  of  the  and  J3 

has  not  been  investigated.  If  the  elements  of  the  perturbing  matrix 
are  small  compared  to  the  elements  of  F  we  may  reasonably  expect 
the  procedure  to  converge  fairly  rapidly. 

§3.5.  Solution  of  Elliptic  Systems  by  the  Big  Matrix: 

Probably  the  most  straightforward  method  for  solving  linear 
elliptic  equations,  provided  that  the  solution  is  specified  on  the 
boundary,  is  to  write  the  complete  set  of  difference  equations  in  the 
form 

<8  Z  =  F  , 

where  Z  is  a  column  vector  of  the  solution,  F  is  a  column  vector 
of  the  boundary  values  and  the  inhomogeneous  term  of  the  P.D.E.  and 
is  the  matrix  of  difference  coefficients.  Following  Bickley  and  McNamee 
[19],  we  will  call  ©  the  big  matrix. 

§5.5.1.  The  big-matrix  representation  of  the  Laplace  and  Poisson 
equations; 

To  illustrate  the  use  of  the  big  matrix  let  us  consider 


Laplace's  equation 
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/afz  A  \  _n 

\  aX2  +  8y2  / 


The  corresponding  difference  equation  is 


-(62  z .  .  +  52  z .  . )  =  c  , 
x  ij  y  ij  ij 


where  c..  is  the  correction  term  and  z.,z  ,  . ,  z,  ,  z,  ,  are 

ij  oj  n+l,j  io  i,m+l 

known. 


Let  us  consider  the  case  of  twelve  interior  points  corres¬ 
ponding  to 

i  =  1,2,3, 
j  =  1,2,3,!+. 

Write  the  difference  equation  in  the  following  way: 


-z.  .  1  -z.  -  .  +4z . .  -z.  .  .  -z.  .  .  =  c..  , 
i,J+l  l-l, j  ij  l+l, J  i,J+l  ij 

i  =  1,2,3;  J  =  1,2, 3, 4. 

This  notation  implies  twelve  equations  for  the  values  of  the  solution 
at  the  twelve  interior  points.  These  equations  can  be  written  in 
matrix  notation 

©Z  =  F  , 

where  Z  is  a  column  of  unknowns*, 

Zt  =  (zirz21>z3l*z12>z22,z32,zl3,Z23,Z33’Zl4’z24,z34^  ’ 

F  is  a  column  of  right-hand  sides, 


Superscript  t  indicates  the  transpose  of  the  matrix. 
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"  ^C11+Z10+Z01,C21+Z20’C31+Z50+Z^1’C12+Z02,C22,C32+Z^2’ 


c13+z03,c13+z03’c23,c33+Z^3’C1^+Z04+z15,c2^+z25, 
c34+z44+Z35}  ’ 


and  (Q  is  the  matrix  of  coefficients, 


Let 


The  big  matrix  can  be  partitioned  into  submatrices  and  written 


“I 


0 


0 

-I 

P3 

-I 


By  following  the  same  procedure  for  general  n 


and  m  we  obtain 
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(Bz  =  F  , 


where 


p 

n 

-I 

0 

0 

...  0 

0 

0 

0 

-I 

P 

n 

-I 

0 

...  0 

0 

0 

0 

0 

0 

0 

0 

...  0 

-I 

p 

n 

-I 

0 

r'O 

0 

0 

...  0 

0 

-I 

p 

n 

Z  =  ^Zll,Z2l,*“',Znl,Z12,Z22,"”,Zn2,""*,Zlm,Z2m,”*,Znm^  » 
Ft  "  (C11+Z10+Z01»C21+Z20 . Cnl+Zn0+Zn+l,l’Cl2+Z02’C22’C52’ 

*  * ' ,Cn-l,2,Cn2+Zn+l,2,Cl3+Z03’ *  *  * ,Cnm+Zn+l,m+Zn,m+l^  * 

©  is  of  order  (nmxnm)  and  Z  and  F  are  of  length  nm.  The 
matrices  P  ,  I  and  0  are  of  order  (nxn).  I  is  the  unit  matrix, 

0  is  the  zero  matrix  and 


P  = 
n 


k  “1  0  ...  0  0  0 

■l  4  »i  ...  o  o  o 


0  0 


0  0 


_0  0  0  ...  0  -1  i 

The  big-matrix  equation  corresponding  to  Poisson's  equation, 


(5.1),  is 


©  Z  =  F  +  G  , 


t  2 

G  =  h  (f11»f21’-”fn1’f12»“*,fnm)  * 


where 
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§5.5-2.  Generalizations  and  discussion; 

Clearly,  any  linear  elliptic  P.D.E.  can  be  written  in  the 
big-matrix  form  and,  provided  a  suitable  linearization  can  be  made, 
quasi-  and  non-linear  elliptic  equations,  too,  are  reducible  to  a 
big-matrix  equation.  Obviously  the  elements  of  the  big  matrix  will 
depend  upon  the  coefficients  of  the  difference  equation,  making  it 
difficult  to  write  down  the  general  form  of  (B  .  The  procedure  for 
setting  up  the  matrix  equation  simply  consists  of  writing  down  the 
difference  equation  at  each  of  the  nm  points  of  the  region  of 
solution . 


It  is  apparent  that  this  method  suffers  from  serious  draw¬ 
backs  which  make  it  of  dubious  practical  value.  Even  with  the  appli¬ 
cation  of  difference-corrections  a  required  accuracy  of  four  or  more 
figures  seldom  makes  consideration  of  fewer  than  one  hundred  interval 
points  practical.  This  means  that  one  hundred  linear  equations  must 
be  solved.  Furthermore,  the  number  of  matrix  elements  to  be  stored 
exceeds  the  memory  capacity  of  most  computers  and,  to  prevent  loss  of 
significance  by  round-off  during  computation,  so  many  figures  must 
be  carried  that  computing  time  would  be  prohibitive.  The  simplicity 
of  the  big-matrix  formulation  masks  its  crudeness.  It  is  a  brutal 
statement  of  the  problem,  but  it  is  the  most  general  closed  difference 
method  and  by  its  crudeness  points  out  how  little  is  known  about 
solving  P.D.  systems. 

§5.6.  Differential  Boundary  Conditions; 

Often  the  boundary  conditions  accompanying  elliptic 
equations  involve  the  normal  derivative  of  the  solution.  Derivative 
conditions  as  such  cannot  be  directly  incorporated  since  we  must  have 


o3  3i<f! i  «  ojtJ.'jpa  ii3*  U  .>  .  .i  on  br»3 
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values  of  the  solution  itself  on  the  boundary.  The  procedure  usually 
followed  is  exemplified  in  the  following  example  involving  Poisson's 
equation*: 


f  dfz  dfz  'N 

W  +  dy2 


f(x,y)  f 


on  the  boundary 


x  =  a  ,  x 
o 


y 


=  b 


1  * 


The  difference  equation  is 


(5.18) 


/  \  ,2, 

-  (5  z.  .  +  &  z  )  =  h  f  .  +  c ,  . 
x  ij  y  ij '  ij  ij 


> 


where  now 


i  =  0»l,2»...,n+l  , 
j  =  0,1,2,.., ,m+l  , 


n  and  m  as  defined  in  §5.1.1,  that  is,  the  difference  equation 
must  be  applied  at  points  on  the  boundary  and  involves  values  of 
the  solution  outside  the  boundary.  The  boundary  conditions  are 


( 


n+1,  j 


dz  N 

57  A  1 

/i,m+l 


=  0  . 


Using  (2.I7), 


h  D  =  p  5  +  c* 

=  \  (E-e'1)  +  c'. 


*  It  is  easily  possible  to  construct  more  complicated  examples 

in  which  considerable  effort  must  be  expended  in  finding  the 
difference  equation  at  nodal  points  on  or  near  the  boundary, 
but  this  would  only  serve  to  obscure  the  main  point  of  this 


section. 
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where  c*  is  a  correction  term,  the  boundary  conditions  are 


z  .  .  =  z ,  +c '  . ,  z  _  .  =  z 
-l»j  ij  OJ  n+2, j 

(5.19) 

,+c'  .  , 

n,j  n+l,j 

ZA  ,  =  Z,,+c','  ,  Z„  _  =  z 

i,-l  il  io  i,m+2 

.  +cV  ,  , 

1m  i,m+l 

where 

c  '  =  2c  *  z 

OJ  X  Oj 

*  Cn+l,J  “  ’ 

2c  '  z  .  , 

x  n+l,j 

c"  =  2c' zs 
io  y  io 

’  Ci,m+1  =  “ 

2c ' z .  . 

y  i,m+l 

We  can  write 

the  difference 

equations  in  matrix  form  in 

a  manner  similar  to  that  of  §5-2.  Let  us  consider  the 

special  case 

n=l 

,  m=2.  The  matrix 

equation  is 

^  1  +  ^2  = 

G  , 

where 

-2-i,o+2zocfzio 

-z-i,i+2zox-zn 

-Z-1,2+2Z02-Z12 

•Z-13+2Z03'Z^ 

ii 

r— 1 

"Z00+2Z10"Z20 

-zoi+2zirz2i 

-Z02+2Z12-Z22 

-z03+2z13'Z23 

'Z10+2Z20'Z30 

•zu+2z2rz3i 

-Z12+2Z22+Z32 

•Z13+2Z23-Z33_ 

~-*o,-i+a*oo-*oi 

-Z00+2Z01-Z02 

"Z01+2Z02'Z03 

■z02+2z03'Z01; 

2  a- 

-zl,-l+2z10'zll 

■Z10+2Zll"Z22 

•Z11+2Z12-Z13 

-z12+2z13-z1H 

_JZ2,-1+2Z20-Z21 

"Z20+2Z2l"Z22 

-Z21+2Z22'Z23 

•Z22+2Z23'Z2lL_. 

G  =  (h2ftj  + 

C^j}  >  i  =  Oj 

1,2;  j  =  0,1,2, 3. 
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The  values  z  . »  z,  , ,  z,  n,  z.  ,  are  replaced  by  the  expressions 
-l,j  3»J  i» "1  v  j  y 

(5.19)  to  give 


2  = 
^1 


2Z00"2Zlo‘C00  2z0l“2zil“c01  2Z02"2Z12“C'o2  2Z03'2z13_C03 
"Z00^2Z10‘Z20  "Z01+2Zll"Z2l  "Z02+2z12~Z22  ~Z03f 2Zl3_Z23 

"2z10+2Z20"C20  '2zll+2z2rC21  “Z12+2Z22"C22  "2z13+2z23_C23 


2-2  0 

-1  2  -1 
0-2  2 


zc 

Z] 

z. 


z  , 

z 

z 

c ' 

c ' . 

c ' 

c ' 

00 

01 

02 

03 

00 

01 

02 

03 

10 

Z11 

Z12 

N 

l—1 

V>l 

- 

0 

0 

0 

0 

20 

Z21 

Z22 

Z23 

C20 

°21 

C22 

C  1 

23 

and 


x2= 


2z00"2z0l"c00  “z00+2z0l“z02  -Z01+2Z02“Z03  ~2Z02+2Z03’C03 

2zl0"2zirCl0  "Z10+2Z11"Z12  “Z11+2Z12'Z13  “2z12+2z13"C13 

2z20"2Z2l"C20  ‘Z20+2Z2rZ22  “Z21+2Z22“Z23  ~2Z22+2Z23~C23 


:oo 

zoi 

Z02 

N 

O 

z 

z 

Z 

:io 

11 

12 

13 

:20 

Z21 

Z22 

Z23 

n 

00 

0 

0 

c"~ 

03 

II 

10 

0 

0 

c" 

13 

II 

'20 

0 

0 

c” 

23 

2 

-2 

0 

0 


-1 

2 

-1 

0 


0 

•1 

2 

■1 


0 

0 

•2 

2 


Let  be  the  matrix  of  order  (kxk)  defined  by 
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D?  = 


~ 2 

-2 

0 

0 

...  0 

0 

~0~ 

-1 

2 

-1 

0 

.  .  O  0 

0 

0 

0 

-1 

2 

-1 

•  •  «  0 

0 

0 

0 

0 

0 

0 

...  -1 

2 

-1 

0 

0 

0 

0 

...  0 

-2 

2 

Then  the  matrix  equation  is 


1)1  Z  +  Z(Dp)t  =  F 


where 


8oo+coo+coo 

g01+Cil 

S02+C02 

803+C03+C 

F  = 

8io+cio 

S11 

812 

Jl20+C20+C20 

g21+C21 

822+C22 

g23+c25+c 

sij  “ h  fu  +  cu  • 


Clearly  the  matrix  equation  for  general  n  and  m  is 
(5.2°)  D?+2  Z  +  Z(D?.,)  =  F  , 


m+2' 


where 


and 


Z  -  U..} 


F  =  fh  fij  +  c!j) 


incremented  in  the  first  row  by  c'.,  in  the  last  row  by  c' 

oj  n+ 1 , j 

in  the  first  column  by  c"  and  in  the  last  column  by  c"  ,,  for 

io  i,m+l 

i  =  0jlj2jooo) n+ 1  j 

j  =  0, 1,2, P . . ,m+l  . 
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Equation  (5-20)  can  be  solved  by  the  methods  of  §5-3*  A 
first  approximation  is  obtained  for  all  correction  terms  equal  to  zero. 
The  solution  is  then  extrapolated  to  give  values  outside  the  boundary, 
correction  terms  are  evaluated  and  included  in  the  right-hand  side  of 
(5-20)  and  the  equation  is  solved  again.  This  procedure  is  continued 
until  successive  solutions  agree  to  the  required  accuracy. 


$5.7.  The  Treatment  of  Singular  Points: 

Up  to  now  we  have  assumed  that  the  boundary  and  initial  data 
are  non-singular .  Three  types  of  singularity  were  discussed  by  Woods 
[21]  in  relation  to  Poisson's  equation,  namely,  discontinuities  in  the 
derivative  of  the  solution  (type  i),  logarithmic  singularities  in  the 
solution  (type  II)  and  simple  discontinuities  in  the  solution  (type 
III).  These  are  probably  the  most  common  types  of  singularity  and 
the  method  of  dealing  with  them  will  exemplify  how  we  may  handle  the 
problem  in  general. 


Suppose  we  are  given  Poisson's  equation, 


c2U  d2U 

ax2  +  ay2 


f(x,y)  > 


with  data  containing  singularities  of  the  types  defined  above.  Since 
the  Laplacian  operator  is  linear  we  can  always  make  the  substitution 


U(x,y)  =  cp(x,y)  +  #(x,y)  . 

The  function  Tp  is  chosen  to  be  a  solution  of  Laplace's  equation  and 

i 

is  such  that  it  contains  the  singularities  leaving  cp  non-singular. 
We  then  solve  the  problem 
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d2cp 

^x2 


+ 


f(x.y) 


with  non-singular  data  in  cp.  The  final  solution,  being  the  sum  of 
<p  and  tp  ,  is  then  calculable  except  possibly,  for  a  neighbourhood 
of  the  singularity. 


In  the  case  of  Poisson's  equation  the  choice  of  tp  is 
relatively  simple  since  we  know  that  the  real  and  imaginary  parts  of 
a  complex  variable  are  solutions  to  Laplace's  equation.  Suppose  the 
point  (r0»^0)  is  a  singular  point.  If  (r0>®0)  is  a  fcyPe  1  singular 
point,  we  choose 

Tp  =  (r-ro)m  (a  sin  m  0  +  3  cos  m  0)m  , 

and,  if  m  is  not  an  integer,  tp  has  a  singular  derivative  of  some 

order  at  (r  ,0  ). 

x  o'  o 

For  a  type  II  singularity  at  (r  ,0  )  we  let 
Tp  =  Re  [log  ({r-r  )e10)m]  =  m  log  (r-r  )  . 

J  O 

Finally,  a  type  III  singularity  determines 
tp  =  Im  [log((r-rQ)  ei0)m]  =  m  0 

A  similar  procedure  is  used  for  other  linear  equations  except 
that  a  function  tp  cannot  always  be  found  which  both  contains  the 
singularities  and  Satisfies  the  homogeneous  equation,  hence,  the  equation 
in  cp  differs  from  the  equation  in  U  in  the  inhomogeneous  term.  This 
is  no  disadvantage  provided  singularities  are  not  introduced  into  the 


equation  by  this  term. 


As  far  as  elliptic  equations  are  concerned,  type  I  and  III 
singularities  cause  no  trouble  if  closed  methods  of  solution  are  used 
and  need  not  be  removed.  Type  II  singularities  must  be  taken  care  of 
in  the  manner  described.  Because  of  the  marching  procedure  used  in 
initial  value  problems  all  varieties  of  singularity  cause  an  instability 
which  is  propagated  in  the  direction  of  increasing  time. 

The  following  example  will  illustrate  the  technique  discussed 

above: 

Suppose  we  are  required  to  solve 

(5»2l)  XUXX  ^xy  +  ^yy  -  ^  * 

U(0,y)  =  | y |  »  (  ^  =  -1,  U(x,+1)  =  1  , 

\  '  x-o 

-1  <  x  <  0  ,  -1  <  y  <  1  . 

The  equation  is  hyperbolic  in  the  region  of  solution  since 

2 

A  =  f  -  x  >  0 
4 

Because 
U(0,y)  =  | y |  , 

there  is  a  type  I  singularity  at  the  origin.  If  we  differentiate 
U(x,+1)  =  1 

with  respect  to  x,  we  have 


whereas  the  second  initial  condition  gives  us 
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y=±i 


Hence,  there  is  a  type  I  singularity  at  the  points 


x  =  0,  y  =  +  1  . 


Let 


u  =  <p  +  , 


where 


=  \  {  *  tan'1  (wi)  •  (if^)  i°s[x2+(i-y2)2]  }  , 


=  |y|  • 


The  system  for  cp  is 


2  f,  r  2  /,  2\21  y^(6y2-3x-8)  +  x(  l+2x)  +  2l 

xcp  -yep  +cp  +“ilog[x  +(l-y  )  ]  +  A  J — 5 - L - - *-z - L -  >- 

XX  'Yxy  yy  It  [  x2  +  ,.2.2  J 


(i-y  )‘ 


<p(0,y)  =  f  (1-y2)  log  (1-y2)  , 

it 


<p(x,+l)  =  -  X  , 

=  -  1  . 

■O 

The  numerical  solution  of  this  system  can  be  undertaken  with  confid¬ 
ence  using  the  method  of  §3.6.2.  The  final  solution  to  (5.21), 


=  cp.  .  + 
ij 


(n>ij 


(*2)iJ 


9 


will  hold  everywhere  in  the  given  region  except  for  a  neighbourhood 


of  the  line  y  =  0  . 
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§5.8.  Non-rectangular  Boundaries: 

If  the  curve  on  which  boundary  data  is  given  is  not  a  rectangle 
with  sides  parallel  to  the  axes, the  methods  discussed  in  the  preceding 
sections  of  this  chapter  are  useless.  Curved  boundaries  do  not  present 
such  a  great  problem  to  relaxation  methods  [3,20]  since  the  mesh  size 
can  be  varied  to  fit  the  boundary.  If  closed  methods  are  used,  the 
usual  procedure  is  to  use  the  given  conditions  to  estimate  data  on  a 
rectangular  boundary.  A  crude  solution  can  thereby  be  obtained  and  a 
better  estimate  of  the  rectangular  conditions  can  be  made  using  this 
solution.  This  process  is  repeated  until  sufficient  accuracy  is  obtain¬ 
ed. 


The  above  technique,  while  proven  by  many  successful  solutions 
to  be  a  valuable  numerical  procedure,  lacks  the  concise  statement  and 
invariant  form  which  distinguishes  closed  methods.  Such  procedures  are 
not  easily  programmed  for  a  digital  computer  since  they  depend,  to  a 
great  extent,  on  the  intuition  and  experience  of  the  computor  rather 
than  a  fixed  set  of  rules.  The  estimation  of  the  first  set  of  condi¬ 
tions  on  the  rectangular  boundary  is  crucial  and  not  easily  undertaken 
by  a  casual  computor. 

In  the  following  section  a  more  straightforward  and  easily 
implemented  method  for  handling  non-rectangular  boundaries  is 
suggested. 

§5.8.1.  The  method  of  linear  combination  of  solutions: 

When  solving  P.D,  systems  analytically  we  often  represent  the 
general  solution  as  a  sum  of  particular  solutions  to  the  equation. 

We  might  reasonably  expect  that  the  solution  to  a  difference  system 
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could  be  obtained  by  combining  a  finite  number  of  solutions  to  the 
difference  equation  in  a  manner  determined  by  the  boundary  conditions. 
The  problem  is  to  obtain  solutions  to  the  difference  equation  and  we 
have  seen  that  this  is  not  possible  in  general  without  referring  to 
boundary  conditions.  Then  if  we  must  solve  a  difference  system,  let 
us  do  so  for  simple  boundary  conditions  on  a  simple  boundary  and 
combine  a  number  of  solutions  obtained  in  this  manner  in  such  a  way 
that  they  satisfy  the  given  conditions  at  a  number  of  points  on  the 
given  boundary.  We  might  then  reasonably  expect  this  combination 
to  be  a  solution  of  the  given  difference  system  and  an  approximate 
solution  to  the  differential  system.  The  most  suitable  boundary 
is  a  rectangle  with  sides  parallel  to  the  axes  and  the  simplest 
boundary  conditions  consist  of  values  of  the  solution. 

In  particular,  let  us  state  the  method  as  follows? 

We  are  given  the  linear  elliptic  equation 

(5.22)  L(Z)  =  0 

and  boundary  conditions 

(5.23)  z  =  zR 

on  the  boundary  B  given  by  g(x,y)  =  0, 

W©  choose  a  rectangular  boundary  R  given  by 

x  *  x  s  a^,  y  s  bQ,  y  b  b^ 

which  encloses  B,  In  R  we  set  up  a  mesh  of  nm  interval  points 
according  to  §5.1.1.  The  difference  equation  is 

( 5 » 2b)  L  '( j )  ”  0  . 
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This  equation  is  to  be  solved  k  times  for  k  different  boundary 
conditions  on  R.  The  choice  of  conditions  is  arbitrary,  but  one 
of  the  more  obvious  is  to  make  the  solution  zero  on  the  boundary 
except  at  one  mesh  point  where  it  takes  on  the  value  unity.  By 
moving  this  point  to  k  different  positions  we  obtain  k  differ¬ 
ent  boundary  conditions.  The  mesh  lines  intersect  R  at  2(m+n) 
points  (corner  points  do  not  appear  in  the  difference  equation  and 
so  are  not  considered).  Hence,  in  this  case 

k  <  2(m  +  n) 

The  matrix  difference  equation  can  now  be  written  down  and  solved 
by  the  methods  of  §5*3»  §5.^  or  §5.5.  We  obtain  k  solutions 

(5.25)  Z(l).  Z(2) .  Z(k)  • 

Let  Z  be  the  matrix  of  the  final  solution.  We  choose 
k  points  on  B  and  call  them 

(5.26)  B^s  » » » »  » 

From  (5.23)  we  have  the  following  known  quantities: 


Each  of  the  solutions  (5.25)  can  be  interpolated  to  the  choosen 


points  (5.26)  on 

B  to  give 

the  k  quantities 

,<0 

V  ’ 

2d) 

»  ...» 

.(i) 

B,  ’ 

1 

2 

k 

z(2) 

z(2> 

.(2) 

B1  • 

>  ...» 

B2 

\  ’ 

,00 

.00 

,(k) 

B1  ’ 

^TJ  9  •  •  0  9 

B2 

Bk  ' 
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We  postulate  constants 


V  a2 . ak 


such  that 


(1)  (2)  (k) 

alZBl  +  *2\  +  •••  +  VBl  "  ZBl  • 


a 


z(l)+  a  z(2)+  .  a  z(k)_ 

1ZB.  2  B„  •'*  +  akZB>  "  B^  * 


(1)  (2)  (k) 

*l\  +  a2\  +  +  ak\  =  \  * 


These  k  linear  equations  in  k  unknowns  can  be  solved  for 


>  ^2 )  •  •  •  f  • 


The  final  solution,  Z,  is  given  by  the  following  linear  combination 

of  the  matrices  (5.25) ‘ 

k 

.(0 


(5.27) 


■  I 


al  z 


i=l 


In  a  sense,  this  is  equivalent  to  the  spectral  resolution 

of  the  matrix  Z  which  we  discussed  in  §5.3.1.  We  let 

n  m 

Z  = 

i=l  j-l 


n  m 

I 1  uij  ri  ^  • 


where  T.  and  cr.  were  known  eigenvectors  and  the  u.  .  were 
i  J  ij 

chosen  to  make  Z  satisfy  the  difference  equation  and  one  set  of 
boundary  conditions.  In  a  similar  manner  we  can  write  (5.27)  in 
the  form 

k  n  m 


z  -  I  I  I  uij<  ri  ai  • 


1=1  i=l  j=l 


100  - 


where  the  u  ^  are  c^osen  to  make  Z  satisfy  the  difference 
equation  and  (k  +  l)  sets  of  boundary  conditions. 

Certain  drawbacks  to  this  method  must  be  pointed  out. 
Computing  time  has  been  increased  considerably  and  storage  for  k 
additional  (nxm)  matrices  is  required.  Some  care  must  be  exercised 
in  the  choice  of  the  conditions  on  R  so  that  we  obtain  independent 
solutions . 


§5.8.2.  Example: 

Consider  Laplace's  equation  in  the  quarter  circle. 


(5.28) 


d2Z  +  dfz 


=  0 


V 


Z(x,0)  =  x\ 

Z(0,y)  =  y\ 

(5.28a)  2  2 

Z(x,y)  =  1  on  x  +  y  =1, 

0  <  x  <  1,  0  <  y  <  1,  x2  +  y2  <  1  . 


The  analytic  solution  is 
(5.29)  Z(r,e)  =  ir*  cos(l+0) 


in  polar  co-ordinates,  or 

Z(x,y)  =  (x2+y2)2  cos(k  tan  1  ^) 

in  Cartesian  co-ordinates.  More  difficult  examples  could  be  chosen 
but  they  would  not  demonstrate  the  numerical  procedure  any  better. 

A  simpler  example  to  investigate  would  have  been  the  system  with 
solution 

r2  cos  (20) 


9 
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but  this  is  also  the  exact  solution  to  the  difference  system  and 
so  would  not  be  a  fair  test  of  the  method. 

Since  the  given  boundary  involves  the  circumference  of  the 
quarter  circle,  application  of  the  method  of  §5.8.1  is  appropriate. 
For  the  rectangular  boundary,  R,  let  us  choose  the  square 

x  =  0,  x  =  1,  y  =  0,  y  =  1. 

2 

Following  the  conventions  of  §5.1.1  we  set  up  a  mesh  of  (n+2) 
boundary  and  interior  points  (x^,  y^),  where 

Xi  =  ih*  yj  = 

i, j  =  0,1,2,... ,n+l, 

1 


h  = 


n+1 


According  to  §5.8.1  we  must  choose  k  sets  of  boundary  conditions 
on  R,  The  matrix  equation  to  be  solved  is  (see  §5.2) 


(5.30)  D2  Z(i)  +  Z^D2  =  F^  , 

'  n  n 


ft  —  1,2,00.,  k 

which  will  give  us  k  solutions.  The  boundary  conditions  on  R 
are  taken  to  be 


(5.30a) 


10 


(ih)U  , 


<*)  „ 

=  (j^r » 
=  8- 


- 

oj 

XO 


i,n+l 


J  i  * 


z«)  .  5 

n+1, j  i , j+n  * 


where 
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o  o  5  2n  *  1  j 

for  i  /  m 


and 

5  =  1  . 

mm 

We  have  taken  k  =  2n°l»  The  figure  below  shows  the  region  of 
solution  for  n  =  b.  Boundary  values  appear  at  their  appropriate 
mesh  points o 


The  first  two  conditions  impose  a  restriction  on  the 

constants  a, »  a  s  „  „  „ ,  a  of  §5<,8»X<>  We  see  this  as  follows 

1  2  2n=>  i 

the  final  solutions,  Zs  is  given  by  (  5°2Y)s> 

2n-l 


ai  Z 


1=1 


But 


.<*> 

•io 


=  zio  =  (ih) 


and 


z(i)  -  *  .  .  (Jh)1* 
oj  oj  'J  ' 
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that  is,  the  final  solution  and  all  the  sub-solutions  have  the  same 
values  along  the  axes.  Hence, 


2n-l 


4=1 


z 


(i) 

io 


=  z 


io 


2n-l 

E 

4=1 


Thus , 


(5-31) 


2n-l 


4=1 


a4  =  1  * 


The  mesh  lines, 

~  ih,  y  =  jh, 
i,j  =  1,2, ... ,n, 

intersect  the  quarter  circle  in  2n  points.  However,  by  choosing 
n  such  that 

(  ,\2  2  .2 

(n  +  1)  =  s  +  t  , 


where  s  and  t  are  integers 
1  <  s,  t  <  n 


we  see  that  there  are  actually  (2n-2)  points  of  intersection. 
Hence,  with  (  5»3l)  we  have  (2n-l)  conditions  on  the  (2n-l) 
constants  a^  .  Labelling  these  points  of  intersection 


B 


2n-2 


(i) 


each  of  the  (2n-l)  solutions,  Z 


9 


can  be  interpolated  to  each 
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2  2 

of  these  points.  Then,  since  the  boundary  values  on  x  +y  =  1 
are  unity,  we  have  that 


2n-l 

I 


J=1 


1  ,  k  =  1,2, . . . , 2n-2, 


2n-l 


i=l 


These  equations  can  be  written  in  matrix  form 
(502)  A  =  B 


and  solved  to  give  the  vector  of  constants 

A  = 

The  final  solution,  (5.27)>  is  then  easily  obtained. 

The  above  procedure  was  programmed  in  Fortran  for  the 

I.B.M.  1620  computer.  Since  the  eigenquantities  are  known  (see 

Appendix  II)  equation  (  5«30)  was  solved  by  the  irrational  method 

of  §5.3.1.  As  sufficient  storage  was  not  available  in  memory  each 

(i) 

sub-solution,  ZK  ,  was  punched  out  on  cards  as  calculated.  Each 
matrix  was  read  back  into  memory  when  called  for  by  the  computer 
and  interpolation  to  the  quarter  circle  was  performed.  Four  terms 
of  a  standard  backward  difference  interpolation  formula  were  used. 
Equation  (5.32)  was  solved  and  the  final  solution  formed.  The 


2n-l 
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numerical  solution  so  obtained  was  compared  with  the  analytic 
solution  (5.29).  The  numerical  solution  and  the  difference  between 
numerical  and  analytic  solution  were  punched  on  cards  and  listed  on 
the  1+07  line  printer.  Dif ference-corrections  were  not  applied. 
Running  time,  as  a  function  of  n,  was  found  empirically  to  be 

t  =  .004n^  +  .In''*  minutes. 

Running  time  for  n  =  9  was  about  one  and  a  half.  The 
solution  for  n  =  9»  to  four  decimal  places,  appears  in  Table 
III.  Values  at  every  second  mesh  point  are  listed.  The  figures 
in  brackets  are  the  differences  between  the  analytic  and  numerical 
solutions . 


TABLE  III 


Approximate  solution  of  the  system  (5.28),  (5.28a). 
Mesh  spacing  =  0.2.  Every  second  value  listed. 


1.0 

1.0000 

(00) 

0.7611+ 

(01) 

0.0651 

(0U) 

0.8 

0.1+096 

(00) 

0.2578 

(-02) 

“0.1786 

(-05) 

-0.8436 

(01+) 

-1.2321 

(02) 

0.6 

0.1296 

too) 

0.0456 

(-08) 

”0.1890 

(-13) 

-0.5174 

(-09) 

“0.81+35 

(03) 

0.4 

0.0256 

too) 

“0.0098 

(-13) 

-0.1007 

(“16) 

“O.1890 

(“13) 

“O.I787 

(-04) 

0.0652 

(03) 

0.2 

0.0016 

(00) 

-0.0053 

(-10) 

”0.0099 

(-12) 

0.01+56 

(-08) 

0.2577 

(-01) 

O.76I3 

(02) 

0.0 

0.0000 

(00) 

0.0016 

(00) 

0.0256 

(00) 

0.1296 

(00) 

0.4096 

(00) 

1.0000 

(00) 

\x 

y  \ 

0.0 

0.2 

0.1+ 

0.6 

0.8 

1.0 
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Considering  the  crudeness  of  the  interval  and  the  fact 
that  no  difference  corrections?  were  used,  the  obtained  accuracy  of 
two  decimals  is  quite  good.  In  order  to  improve  significantly  the 
accuracy  of  the  solution,  difference-corrections  must  be  applied. 
Because  of  the  fourth  order  dependence  of  time  on  mesh  size,  re¬ 
duction  of  mesh  spacing  to  improve  accuracy  would  be  prohibitively 
costly  in  time. 

It  is  clear  that  this  method  for  handling  non-rectangular 
boundaries  requires  further  investigation.  Rates  of  convergence  and 
the  effect  of  the  choice  of  matching  points  on  the  given  boundary 
should  be  considered.  It  is  likely,  however,  that  this  procedure 
can  be  used  with  about  as  much  confidence  as  any  other  closed  diff¬ 
erence  method.  Boundary  conditions  involving  derivatives  of  the 
solution  can  be  similarly  handled. 


.... 
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CHAPTER  VI 

CONCLUSION 

In  the  preceding  chapters  an  attempt  has  been  made  to  describe 
how  one  goes  about  solving  a  second  order  partial  differential  system 
in  two  independent  and  one  dependent  variables.  Much  of  what  has  been 
said  can  be  applied  to  equations  of  higher  order  and  more  independent 
variables.  In  the  way  of  a  general  method,  all  that  can  be  said  is  that 
the  differential  system  is  changed  to  a  difference  system  by  replacing 
derivatives  by  terms  from  a  difference  series  and  the  difference  system 
is  solved  numerically.  This  procedure  has  been  illustrated  briefly  in 
reference  to  hyperbolic  and  parabolic  equations  (Chapters  III  and  IV) 
and  in  more  detail  with  respect  to  elliptic  systems  (Chapter  V). 

i 

Examples  have  been  given. 

Except  for  a  few  special  cases,  the  numerical  solution  cannot 
be  justified  by  a  priori  arguments.  We  expect  that  trouble  will  mani¬ 
fest  itself  in  some  obvious  way  such  as  instability  or  non-convergence 
of  successive  approximations  to  the  solution  (Chapter  IV). 

Matrix  methods  for  solving  elliptic  systems  have  been  discussed 
in  detail  (Chapter  V).  Where  applicable,  these  methods  provide  a  stra¬ 
ightforward  and  easily  implemented  algorithm  for  obtaining  a  solution. 

On  the  whole,  each  partial  differential  system  must  be  consid¬ 
ered  on  its  own  merits  which  clearly  indicates  that  our  knowledge  of 
partial  differential  equations  -  analytical  and  numerical  -  is  far  from 


complete. 
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